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Abstract 


The  finite  element  method  was  used  to  solve  the  nonlinear, 
small -disturbance,  transonic,  velocity-potential  equation  for  problems 
of  steady  flow  over  a  circular  cylinder  and  over  a  thin-airfoil  in  a 
uniform  steady  airstream.  The  governing  differential  equation  is  valid 
for  inviscid,  irrotational ,  isentropic  flow  of  a  perfect  gas  to  in¬ 
clude  weak  shocks  providing  airflow  separation  does  not  occur.  For 
compressible  subsonic  and  transonic  flows  the  nonlinear  small -disturbance 
equation  was  expressed  in  iterative  form  as  a  sequence  of  linear  equa¬ 
tions  which  was  solved  iteratively  until  the  difference  between  two 
successive  solutions  became  arbitrarily  small.  For  analysis  purposes 
the  infinite  flowfield  was  replaced  by  a  finite  but  sufficiently  large 
domain  that  was  discretized  with  sector  elements  for  the  cylinder  pro¬ 
blem  and  rectangular  elements  for  the  airfoil  problem.  The  finite 
element  equations  were  obtained  from  Galerkin's  Method  of  Weighted 
Residuals.  Boundary  conditions  of  the  Neumann  type  were  imposed  along 
the  surface  contour  of  the  cylinder  and  along  an  approximate  boundary 
in  accordance  with  classical  thin-airfoil  theory  for  the  airfoil.  For 
both  problems  Dirichlet  conditions  were  imposed  along  the  farfield 
boundary  from  an  asymptotic  solution  which  satisfies  the  actual  infinity 
condition  and  is  valid  in  the  farfield. 

Three  sub-problems  were  investigated  for  the  circular  cylinder. 
First,  three  different  types  of  trial  functions  were  investigated  to 
approximate  the  solution  for  the  velocity  potential  function  for  the 
case  of  incompressible  flow  without  circulation.  The  three  trial 
functions  were:  (1)  a  trignometric  approximation  resulting  in  a 


x 


non-conforming  element,  (2)  a  bilinear  polynomial  (conforming  element) 
typical  of  elements  used  in  finite  element  analyses,  and  (3)  a 
rational  approximation  resulting  in  a  new  conforming  element.  Con¬ 
vergence  properties  of  each  element  were  studied  as  a  function  of  dis¬ 
cretization  refinement  (element  size) .  The  new  element  proved  to  have 
superior  properties  for  the  problem  solved. 

The  second  subproblem  for  the  cylinder  was  to  use  the  two  conform¬ 
ing  elements  to  obtain  solutions  for  incompressible  flow  with  circulation. 
Superposition  was  used  to  split  the  total  problem  into  two  elementary 
component  problems.  The  value  of  circulation  was  determined  by  en¬ 
forcing  the  stagnation  or  Kutta  condition  after  each  component  solution 
was  found.  An  optimum  way  to  select  and  refine  the  discretization  was 
discovered  to  insure  that  the  error  in  circulation  was  kept  to  a 
minimum. 

The  third  subproblem  for  the  cylinder  was  to  solve  the  compressi¬ 
ble  flow  problem  without  circulation  using  the  new  element.  In  the 
strict  sense,  the  small-disturbance  equation  is  not  valid  for  compressi¬ 
ble  potential  flow  over  a  cylinder;  however,  the  problem  is  mathematically 
well-posed  with  its  use.  The  iterative  solution  scheme  converged 
rapidly  for  subsonic  flows,  but  failed  to  converge  for  transonic  flows 
when  the  supersonic  zone  engulfed  one  or  more  complete  elements.  Pre¬ 
dictions  of  the  critical  Mach  number  and  subsonic  velocity  distributions 
are  compared  with  known  results. 

Convergence  of  iterative  finite  element  solution  schemes  for 
transonic  (mixed-elliptic-hyperbolic)  flows  was  investigated  for  the 
airfoil  problem.  An  examination  was  made  of  the  divergence  behavior 
of  solution  schemes  reported  by  other  investigators  who  have  used 


xi 


finite  element  formulations  similar  to  those  used  in  this  study.  As  a 
result  of  this  examination  an  iterative  solution  procedure  was  developed 
which  permitted  convergence  of  solution  schemes  for  mixed  (transonic) 
flows.  The  procedure  includes  a  new  "upwinding"  technique  that  accounts 
for  the  proper  zone  of  influence  for  elements  in  the  supersonic 
(hyperbolic)  region.  Governed  by  two  parameters,  the  new  technique 
alters  the  finite  element  formulation  to  exclude  the  influence  of  itera¬ 
tive  downwind  forces  on  the  solution  at  upwind  nodes.  The  "upwinding" 
technique  not  only  "arrests"  the  divergence  behavior  of  the  solution 
scheme,  but  also  "captures"  the  weak  compression  shock  which  forms 
automatically  without  the  use  of  shock  elements. 

Known  results  from  experimental  data,  classical  solutions,  finite 
difference  solutions,  and  other  finite  element  solutions  are  compared 
with  the  finite  element  solutions  obtained  in  this  study.  Comparisons 
of  velocity  and  pressure  distributions  are  given  for  cases  of  incompressi 
ble  flow  and  compressible  subsonic  and  transonic  flows. 


PREDICTION  OF  AERODYNAMIC  FORCES  ON  A  CIRCULAR  CYLINDER 


AND  A  THIN  AIRFOIL  IN  A  TRANSONIC  AIRSTREAM 
BY  THE  FINITE  ELEMENT  METHOD 


I  Introduction 


Analysis  of  transonic  flow  is  one  of  the  most  challenging  problems 
in  potential  aerodynamics  today  due  to  the  nonlinear  nature  and  mixed 
character  of  the  flow.  Transonic  flow  may  occur  with  modern  aircraft 
during  flight  maneuvers,  encounters  with  atmospheric  turbulence  or  wind 
gusts,  and  during  accelerations  to  supersonic  speeds.  Many  fighter 
aircraft,  for  example,  often  endure  extensive  portions  of  their  mission 
profile  at  transonic  speeds.  During  these  periods  violent  oscillatory 
motion  may  occur  which  could  pose  a  hazard  to  flight.  Thus,  the  need 
is  evident  for  accurate  and  reliable  methods  of  analysis  to  predict 
aerodynamic  loads  at  transonic  speeds. 

A  renewed  interest  in  transonic  flow  is  manifest  by  the  rather 
large  volume  of  technical  papers  that  have  appeared  in  the  literature 
within  the  last  fifteen  years.  Several  methods  of  analysis  have  been 
developed,  but  perhaps  those  most  extensively  used  today  are  finite 
difference  methods.  The  majority  of  these  methods  are  devoted  to  two- 
dimensional  potential  flow  analysis,  although  the  flow  is  more  complex 
in  nature.  This  simplifying  position  is  justifiable  to  some  degree 
since  selection  of  suitable  section  shapes  has  always  been  one  of  the 
stages  in  the  design  process  for  aircraft  wings  and  helicopter  blades. 
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For  design  purposes,  and  for  many  other  problems  of  engineering 
interest,  the  flow  about  a  body  may  be  adequately  described  by  potential 
aerodynamics  and  small  disturbance  theory  (Ref  1).  Consequently,  a 
large  portion  of  the  technical  papers  on  transonic  flow  adopt  this 
simplifying  position.  The  assumptions  of  two-dimensional,  inviscid, 
irrotational ,  isentropic  flow  drastically  simplify  the  coupled  basic 
equations  of  fluid  mechanics.  Introducing  a  velocity  potential  produces 
a  single,  governing,  second-order,  nonlinear,  partial  differential 
equation  which  is  valid  for  flows  with  weak  shocks.  Small  disturbance 
approximations  provide  further  reduction  in  analytical  complexity. 

Errors  in  pressure  distributions  resulting  from  these  assumptions  are 
not  severe  except  at  stagnation  points,  at  large  angles  of  attack,  or 
when  extensive  airflow  separation  occurs. 

Within  the  last  eight  years  a  relatively  new  method  of  analysis  has 
been  vised  to  solve  potential  aerodynamics  problems .  This  new  method 
was  originally  developed  by  structural  engineers  and  is  known  as  the 
Finite  Element  Method  (FEM) .  Its  application  to  airfoil  analysis  is 
currently  being  investigated  by  both  the  aircraft  industry  and  govern¬ 
mental  agencies. 

Previous  and  Recent  Works 

Perhaps  the  first  paper  to  propose  the  use  of  finite  elements  for 
field  problems  involving  Laplace  and  Poisson  equations  was  by 
Zienkiewiez  and  Cheung  in  1965  (Ref  2) .  It  was  three  years  later  before 
the  FEM  was  used  to  solve  an  aerodynamics  problem.  Martin  (Ref  3)  used 
linear  triangular  elements  and  a  variational  principle  to  solve  for  the 
stream  function  for  incompressible  flow  over  a  circular  cylinder  located 


between  parallel  walls.  Within  a  year  Norrie  and  deVries  (Ref  4-8) 
developed  finite  element  techniques  to  solve  incompressible  problems 
for  flow  over  single  and  cascading  airfoils.  They  formulated  the 
problem  in  terms  of  the  velocity  potential  and  used  linear  triangular 
elements  along  with  a  variational  principle  to  obtain  the  finite  element 
equations.  They  also  made  liberal  uses  of  the  superposition  principle. 
For  example,  the  problem  of  lifting  flow  over  an  airfoil  was  divided 
into  thickness  and  lifting  parts,  each  of  which  was  solved  by  the  FEN. 

The  solutions  were  linearly  combined  and  the  Kutta  condition  applied 
to  specify  the  circulation.  Unfortunately,  they  were  not  able  to  present 
many  computational  examples  due  to  computer  limitations. 

Shen  (Ref  9)  published  an  interesting  paper  with  intent  to  "bring 
the  maximum  amount  of  a  priori  information  theoretical  and  analytical, 
so  as  to  minimize  the  chore  that  must  be  done  numerically  in  the  Finite 
Element  Method."  He  formulated  the  problem  of  incompressible  flow  over 
a  lifting  airfoil  in  terms  of  the  stream  function  and  vised  a  variational 
principle.  Shen  modeled  the  infinite  domain  with  an  inner  and  outer 
patching  technique  similar  to  that  used  in  finite  difference  relaxation 
schemes  (Ref  10) .  The  infinite  domain  was  divided  into  two  super¬ 
elements  or  patches.  The  inner  patch  contained  the  airfoil  and  a  portion 
of  the  flow  field  extending  from  the  airfoil  to  some  arbitrary  but 
sufficiently  large  distance  from  it.  The  outer  patch  completed  the 
infinite  flowfield  within  which  an  analytically  obtained  asymptotic 
solution  with  undetermined  parameters  was  used.  The  FEN  was  used  to 
obtain  the  solution  only  in  the  inner  patch.  Globally  the  two 
solutions  were  matched  along  the  comnon  boundary  separating  the  twa 
patches.  Shen  presented  some  results  for  the  circular  cylinder  and 


Joukowski  airfoils  using  linear  triangular  elements,  but  stated,  "in 
actual  implementation  of  the  finite  element  method  for  the  inner  patch, 


many  details  remain  that  affect  the  accuracy  of  the  results." 

For  steady,  compressible,  potential  flow  over  an  arbitrary  body 
additional  difficulties  appear.  The  exact  governing  equation  is  non¬ 
linear,  and  when  the  free  stream  Mach  number  becomes  large  enough, 
bubbles  of  supersonic  fiow  may  appear  over  a  portion  of  the  body.  For 
totally  subsonic  flow  Periaux  (Ref  11)  was  successful  in  using  linear 
and  quadratic  triangular  elements  combined  with  an  iterative  solution 
algorithm  to  solve  for  potential  flow  over  airfoils.  Solutions  were 
obtained  from  both  velocity  potential  and  stream  function  formulations 
by  minimization  of  functionals.  Shen  and  Habashi  (Ref  12)  noted  that 
solutions  derived  from  variational  principles  would  not  converge  for 
supercritical  Pfech  numbers.  They  proposed  a  local  linearization  of 
the  problem  which  could  be  formulated  in  terms  of  either  the  stream 
function  or  the  velocity  potential  function.  The  local  governing  dif¬ 
ferential  equation  for  the  local  perturbation  velocity  potential  $ 
was  shown  to  be  the  linear  small-disturbance  equation  ( I  ~  Mg  j  $ff 
4»  rh  —  rt  .  The  elemental  Mach  number  M  was  taken  to  be  con- 

stant  in  element  £  and  was  calculated  from  the  previous  iteration. 
Coordinates  (M)  were  aligned  parallel  and  normal  to  the  previous 
iteration  of  the  streamlines.  They  then  used  linear  triangular  elements 
and  the  Prandtl-Glauert  transformation  to  solve  for  compressible  flow 
over  a  circular  cylinder  without  circulation.  Flcwfields  for  airfoil 
problems  were  discretized  by  a  mapping  procedure  using  the  inverse 
Joukowski  transformation.  Their  results  converged  for  compressible 
flow  to  include  a  small  supersonic  bubble,  but  were  valid  only  for 


Mach  numbers  near  the  critical  values. 

To  date  the  most  extensive  application  of  the  FEM  for  airfoil 
analysis  was  done  by  Chan  and  Brashears  (Refs  13-17) .  In  their  initial 
work  the  small-disturbance  velocity  potential  equation  was  solved  for 
the  cases  of:  incompressible  flow  about  a  lifting  airfoil,  steady 
compressible  subsonic  and  transonic  flows  about  a  nonlifting  airfoil, 
and  unsteady  transonic  flow  about  an  airfoil  that  is  harmonically 
oscillating  about  a  nonlifting  mean  state.  The  finite  element  solution 
was  obtained  with  cubic  triangular  elements  using  Galerkin's  method 
of  weighted  residuals.  For  steady  compressible  flow  the  governing 
differential  equation  was  cast  into  an  equation  of  the  Poisson  type  for 
which  finite  element  equations  were  constructed  and  solved  by  iterative 
algorithms.  The  infinite  domain  was  divided  into  two  sub-domains  or 
patches  as  described  previously.  A  solution  was  obtained  in  the  inner 
patch  by  finite  element  techniques  with  farfield  boundary  conditions 
specified  by  the  farfield  expressions  of  Klunker  (Ref  18) .  The  un¬ 
steady  transonic  flow  problem  was  treated  as  a  sum  of  two  problems  after 
Landahl  (Ref  19) :  (1)  the  non-linear  transonic  flow  problem  about  a 

mean  steady  position,  and  (2)  a  linear  oscillation  problem  about  the 
mean  steady  position.  Generally  speaking,  velocity  and  pressure  distri¬ 
butions  compared  well  with  other  analytical  or  experimental  data  for 
nonlifting  subsonic  flows.  For  steady  transonic  flow,  the  iterative 
Galerkin  formulation  failed  to  converge  and  was  abandoned  in  later 
efforts.  For  lifting  cases  the  finite  element  solutions  did  not 
compare  as  well  with  known  results  as  solutions  did  for  nonlifting  flows. 
The  major  difficulties  encountered  were  determining  correct  values  of 
circulation  and  accurate  estimates  of  the  pressure  distributions  near 
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the  leading  and  trailing  edges.  Predictions  of  pressure  distributions 
for  unsteady  flows  were  only  quantitatively  comparable  with  experimental 
results . 

In  later  works  Chan  and  Brashears  (Refs  15-17)  abandoned  the 
Galerkin  formulation  used  in  their  earlier  work.  They  tried  several 
techniques  to  eliminate  the  improper  downwind  influences  upon  the  solu¬ 
tion  at  upwind  supersonic  nodes,  but  were  unsuccessful  when  the  Galerkin 
method  was  used.  They,  instead,  adopted  the  least  squares  method  of 
weighted  residuals,  but  unfortunately  used  elements  which  were  not 
compatible  with  the  new  formulation. 

Recently,  alternative  implicit  velocity  formulations  for  small- 
disturbance  theory  have  been  suggested  by  Wellford  &  Hafez  (Refs  20-21) , 
Ecer  (Ref  22) ,  and  Aziz  (Ref  23) .  In  order  to  make  solutions  converge , 
Wellford  suggested  adding  time  dependent  and  explicit  artificial  vis¬ 
cosity  terms  to  the  governing  differential  equations,  expressed  in 
terms  of  the  velocity  perturbations.  Akay  (Ref  24)  presents  a  finite 
element  model  for  the  full,  two-dimensional,  potential  equations  using 
variational  principles.  Solutions  for  subcritical  flows  showed  good 
agreement  with  experimental  data,  but  failed  for  mixed  flows.  The  use 
of  artificial  viscosity,  both  explicit  and  implicit,  was  introduced  in 
the  analysis  to  prevent  divergence  of  solutions.  The  addition  of 
viscosity  caused  an  oscillation  about  some  solution  in  the  solution 
algorithm  used,  but  it  did  not  provide  convergence  in  the  absolute  sense. 

Obj  ective 


The  purpose  of  this  work  was  to  use  the  finite  element  method 
(FEW)  to  predict  surface  pressure  distributions  for  two-dimensional, 


potential  flow  over  bodies  in  an  infinite  uniform  flowfield  for  incom¬ 
pressible  and  compressible  flows  that  include  the  transonic  regime  with 
weak  shocks.  The  Galerkin  method  of  weighted  residuals,  abandoned  by 
Chan  and  Brashears,  was  used  with  appropriate  conforming  "lower"  order 
elements.  The  principle  objective  of  this  study  was  to  make  as  many 
simplifying  assumptions  as  possible  in  both  the  flow  model  and  the 
numerical  approximation  procedure,  and  determine  whether  a  "simple" 
Galerkin  approach  is  acceptable  for  solving  a  rather  complicated,  non¬ 
linear,  mixed-flow  problem  in  an  uiiuounded  region  where  discontinuities 
and  singularities  may  exist. 

This  study  was  divided  hlto  two  parts,  flow  over  a  circular 
cylinder  and  flow  over  i:  \hin  airfoil.  One  of  the  original  intents  of 
the  stuay  was  to  examine  some  of  the  characteristics  of  the  FEM's 
application  to  potential  flow  problems.  Although  this  intent  was 
modified,  it  was  accomplished,  in  part,  with  the  circular  cylinder 
problem.  Three  different  elements  were  used,  representing  conforming 
and  non-conforming  elements.  Trial  solutions  for  bilinear,  rational, 
and  trignometric  approximations  were  examined.  Properties  of  flow 
field  discretization  and  application  of  boundary  conditions  were 
examined. 

The  numerical  procedure  for  solving  the  velocity  potential 
equation  for  transonic  flows  was  confined  to  the  airfoil  problem. 

Small  disturbance  theory  and  approximations  from  classical  thin- 
airfoil  theory  were  used.  The  purpose  for  examining  the  transonic 
regime  was  to  investigate  convergence  problems  reported  by  other 
investigators,  and  attempt  to  develop  a  technique  to  account  for  the 
proper  zone  of  influence  in  the  supersonic  region.  If  an  adequate 
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method  could  be  developed  to  exclude  the  downwind  influence  on  the 
solution  at  upwind  supersonic  nodes,  then  the  Galerkin  procedure  would 
not  have  to  be  abandoned  as  concluded  by  Chan  and  Brashears. 

The  Finite  Element  Method  (Background  Information) .  The  finite  element 
method  was  developed  about  25  years  ago  by  structural  engineers  to 
analyze  complex  structural  systems.  An  engineering  structure  can  be 
thought  of  as  a  collection  of  discrete  elements  that  are  inter-connected 
at  a  finite  number  of  locations.  For  example,  a  simple  pin-connected 
truss  can  be  modeled  by  axial-rod  elements.  For  a  given  static  loading, 
the  enforcement  of  global  equilibrium  (force  balance)  is  sufficient  to 
determine  the  tensile  or  compressive  loading  of  each  element.  In  a 
continuum,  whether  it  be  structural  or  otherwise,  the  number  of  con¬ 
nections  becomes  infinite.  For  such  problems  the  continuum  must  be 
divided  into  a  finite  number  of  conveniently  shaped  elements  and  trans¬ 
formed  into  a  discretized  finite  assemblage  of  nodal  parameters. 

In  the  earlier  stages  of  finite  element  development  static  force 
balance  formed  the  theoretical  foundation  of  the  method.  Later,  energy 
principles,  which  form  a  significant  part  of  structural  analysis,  were 
used.  The  finite  element  method  provided  a  way  to  approximate  the 
global  strain  energy  of  a  continuous  structural  system  in  terms  of  the 
ensemble  of  energies  in  local,  discretized  subsystems  or  elements.  Nodal 
displacements  are  determined  frcm  admissible  assumed  displacement  dis¬ 
tributions  by  minimizing  the  strain  energy  functional. 

In  the  mid-1960s  the  FEM  was  examined  for  possible  application  to 
non-structural  problems,  such  as  fluid  flow.  Successful  application 
led  to  further  development  of  the  theory.  It  has  since  been  generalized 
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to  solve  linear  and  nonlinear  partial  differential  equations  for 
boundary  and  initial  value  problems  in  many  fields  of  engineering  and 
mathematical  science.  Recently,  several  mathematicians  have  contri¬ 
buted  to  the  Finite  Element  development  and  have  established  it  as  an 
important  branch  of  approximation  theory.  Close  relationships  exist 
between  finite  element  analysis  and  the  classical  variational  concept 
of  the  Rayleigh-Ritz  method  in  problems  where  variational  principles 
apply.  Unfortunately,  variational  principles  cannot  be  found  in  all 
engineering  problems,  particularly  when  governing  differential  equations 
are  not  self-adjoint.  However,  for  such  problems  the  weighted  residual 
methods,  such  as  least  squares,  collocation,  or  the  well-known  Galerkin 
method  apply.  The  Galerkin  method  is  perhaps  the  most  convenient 
weighted-residual  method  for  FEM  analysis. 
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II  Problem  Descriptions 

Two  steady,  potential  flow  problems  are  considered  in  this  study. 
The  first  is  flow  about  a  circular  cylinder  and  the  second  is  flow  over 
a  thin  symnetric  airfoil  at  zero  angle  of  attack. 


Circular  Cylinder/Incompressible  Flow 

First  consider  the  problem  of  steady,  incompressible,  inviscid, 
irrotational  flow  about  a  unit  circular  cylinder  (i.e.  radius  =  1) 
placed  in  a  uniform  steady  airstream  of  infinite  extent.  Let  the  free- 
stream  be  directed  in  the  positive  X -direction  with  coordinate  systems 
attached  to  the  center  of  the  cylinder  as  shown  in  Fig  1.  Let  SL 
denote  the  infinite  flowfield  domain  composed  of  points  (f,  6)  .  The 
boundary  of  Si  (denoted  d-H.  )  is  composed  of  all  points  ( f,  & )  on  the 
cylinder  surface  (denoted  )  and  the  boundary  at  infinity  (denoted 

d&ep  ,  i.e.  points  (fJS.)  as  - - -  CO  )  .  For  flow  with  circulation 

a  branch  cut  is  placed  in  XI  as  shown  in  Fig  1.  The  purpose  of  the 
cut  is  discussed  in  a  subsequent  section. 

Governing  Equation.  The  governing  differential  equation  for  in¬ 
compressible  potential  flow  is  the  Laplace  equation  given  by 

O  for  C  f,6)i-n  fl  (1] 

where  £  is  the  velocity  potential  function.  Since  the  equation  is 
linear,  the  potential  function  can  be  expressed  as 
The  term  the  potential  of  the  free -stream  with  velocity  \feo 

and  the  perturbation  velocity  potential.  Substituting  this 
expression  into  eq  1  gives 


(2) 


mm 


^z<f>  (r,e)-  o  for  (^0)  in  -ft 

which  is  the  governing  equation  for  (j)  . 

Boundary  Conditions .  For  the  problem  to  be  well  posed,  boundary 
conditions  must  be  specified  on  the  surface  of  the  cylinder  and  at 
infinity.  At  infinity  the  disturbance  velocities  must  vanish 

v  »  O  as  r— (3) 

The  tangential  condition  applied  along  the  surface  of  the  cylinder  is 

y  §  •  H  -  O  for  (r,  9)  in  ^ Sic. 

where  A  is  a  unit  vector  pointing  outward  from  the  surface.  This 
condition  reduces  to 

(for  +  C .OS  6  —  C>  for  (r,  $)  in  (5) 

For  flow  with  circulation  the  stagnation  condition  must  be 
enforced  at  the  down-stream  stagnation  point,  which  corresponds  to 
enforcing  the  Kutta  condition  at  the  trailing  edge  of  an  airfoil. 

This  condition  is  given  by 

-  o  at  (r,e)-  (i,-?*)  (6) 

Since  -  O  from  the  tangential  condition,  then  the  stagnation  or 
Kutta  condition  is  satisfied  by 

(<pie-  rsine)  |  -  ° 
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In  order  to  keep  <p  single  valued,  a  branch  cut  is  placed  in  the  flow- 
field  as  previously  mentioned.  It  can  be  shown  that  the  jump  in 
potential  across  the  cut  ( r)  is  equal  to  the  circulation 
defined  by 

F  =  §  v  •  d  s  (8) 

V  is  the  tangential  velocity  and  ds  is  tangent  to  the  cylinder.  Thus, 
across  the  cut  the  condition 

r  -  </>*-  <f>~  (9) 

must  be  enforced.  The  actual  value  of  P  is  unknown,  but  can  be  deter¬ 
mined  by  enforcing  the  stagnation  condition.  The  numerical  procedure 
for  doing  this  is  described  in  the  next  chapter. 

When  the  entire  domain  JTL  is  not  discretized,  then  an  additional 
boundary  condition  is  needed  when  symmetric  flow  (no  circulation)  is 
considered.  The  problem  of  flow  without  circulation  can  be  sol  -td  ;.n 
the  upper -half  space  defined  by  y  »  O  (denoted  -ft /a) .  This  is  possible 
due  to  symmetry.  Actually,  any  quadrant  of  H  would  be  sufficient  to 
solve  for  ^  in  all  of  SL  ■  In  the  half-space  the  additional  boundary 
condition  imposed  along  the  axis  y~o  for  all  points  f  >  I  (denoted 
£1*0  )  is 

~  &  along  (10) 

This  condition  is  enforced  from  symmetry  considerations  and  is  not  a 
physical  boundary  condition. 


13 


Thin  Airfoil 


Consider  the  problem  of  steady,  inviscid,  irrotational ,  isentropic 
flow  about  a  thin  airfoil  placed  in  a  steady,  uniform  airstream  of  in¬ 
finite  extent .  Let  the  free  stream  be  directed  in  the  positive  X  - 
direction  with  the  airfoil  chord  aligned  with  the  X-axis.  Let  XL 
denote  the  infinite  flowfield  dcmain  composed  of  points  OM)  • 

The  boundary  of  XL  (denoted  dXL  )  is  composed  of  all  points  (X,  y  ) 
on  the  airfoil  surface  (denoted  dXL^)  and  the  boundary  at  infinity 
(denoted  i.e.  points  )  as  f — >00),  as  depicted  in 

Fig  2.  For  cases  involving  lift  a  cut  is  placed  in  the  flowfield 
leading  from  the  trailing  edge  to  the  boundary  at  infinity.  The 
reasons  for  the  cut  are  essentially  the  same  for  the  airfoil  as  for 
the  cylinder  with  circulation,  which  has  been  previously  described. 

Governing  Differential  Equation.  For  a  thin  airfoil,  small- 
perturbation  theory  can  be  used  to  describe  the  disturbances  in  the 
free-stream  velocity  caused  by  the  presence  of  the  airfoil.  The 
governing  differential  equation  for  such  a  problem  is  the  well-known 
small-disturbance  equation  for  the  non-dimensional  velocity  potential 
function  (Ref  19:4).  The  equation  is 

[|-mL-mL0+  *)  A xx  i  iyy  r  0  (11) 

for  all  points  (%,y)  in  fl  .  Mae  is  the  free  stream  Mach  number  and 
IS  the  ratio  of  specific  heats  ("&  =1.4  for  air).  The  dimensionless 
variables  are  related  to  the  physical  ones  (  X,7) 

by  the  relations :  X  '  *  /  C  j  y  ~  V/d  ;  4>  ~  $> /  \Joo  C. 


Figure  2  -  Flowfield  for  the  Airfoil 


where  C  is  the  airfoil  chord  length  and  \Tce  »  the  free  stream  velocity. 
When  considering  strictly  transonic  flow  over  an  airfoil  with  a  given 
thickness  ratio,  slightly  different  transformations  relating  y  and  y 
are  vised  (Ref  10) .  In  the  strict  sense ,  eq  11  is  valid  only  for  isen- 
tropic  flow,  but  it  is  a  good  approximation  for  flows  with  weak  shocks 
(Ref  19:2). 

For  incompressible  flow  (M„  =  0)  ,  eq  11  reduces  to  the  Laplace 
equation  which  is  valid  for  bodies  of  arbitrary  shape  (e.g.  the  circular 
cylinder).  For  low  subsonic  flow  (i.e.  the  nonlinear 

term  °ften  assumed  small  compared  to  the  other  terms  and 

dropped,  as  done  in  linear  theory.  However,  for  flows  in  the  transonic 
range,  the  nonlinear  term  becomes  important  and  must  be  kept  to  model 
the  mixed  subsonic-supersonic  flow  regime  accurately.  In  this  study 
the  nonlinear  term  is  retained  even  for  low  subsonic  flow. 

Equation  11  is  a  nonlinear,  second-order,  partial  differential 
equation  of  mixed  type.  The  mixed  character  is  caused  by  a  change  in 
sign  of  the  coefficient.  For  subsonic,  free-stream  Mach  numbers 

the  sign  change  is  caused  by  the  behavior  of  the  nonlinear  term.  For 
subcritical  values  of  the  variable  coefficient  of  <$,xX  is  positive 
and  eq  11  is  elliptic  everywhere  in  £1  .  When  is  increased 

slightly  beyond  some  critical  value,  the  coefficient  becomes  negative 
in  a  region  in  the  vicinity  of  the  maximum  thickness  point  of  the 
airfoil.  Equation  11  is  hyperbolic  in  this  region  while  it  remains 
elliptic  in  the  remainder  of  JCZ,  .  Along  the  line  which  separates  these 
mixed  flow  regions,  the  coefficient  is  zero  and  the  equation  is  para¬ 
bolic.  Unfortunately,  the  position  of  the  parabolic  or  sonic  line  is 
not  known  a  priori.  In  addition,  the  existence  of  the  nonlinear 


term  permits  solutions  which  have  discontinuous  first  derivatives. 

These  discontinuities  are  associated  with  the  presence  of  weak  com¬ 
pression  shocks  which  separate  the  downstream  side  of  the  hyperbolic 
and  elliptic  regions.  These  shocks  permit  the  supersonic  flow  in  the 
hyperbolic  region  to  return  to  subsonic  flow  in  the  elliptic  region 
over  very  small  distances. 

Boundary  Conditions .  For  the  problem  to  be  well  posed,  a  boundary 
condition  must  be  specified  on  the  surface  of  the  airfoil  and  at 
infinity.  At  infinity  the  disturbance  velocity  must  vanish,  i.e. 


o  as  r  - (12) 

The  impervious  or  tangential  condition  for  flow  along  the  airfoil  sur¬ 
face  is  given  by 

=  ( i  +  h)  for  ^  <13) 

where  V-'fC*)2  O  describes  the  surface  of  the  airfoil. 

To  be  consistent  with  linear  small-perturbation  theory,  the  4  X  term 
in  the  tangential  boundary  condition  is  normally  neglected  compared  to 
unity.  This  boundary  condition  is  not  enforced  on  the  surface 

Instead,  it  is  applied  along  the  axis  yz  0~  >  in 
accordance  with  classical  thin-airfoil  theory  (Ref  58) . 

For  nonlifting  thin  airfoils  the  solution  can  be  obtained  in  the 
upper  half-space.  From  synmetry  considerations  the  condition 

( fa  =  O  (14) 

must  be  applied  for  all  points  along  the  axis  >fsO  which  lie  beyond 
the  leading  and  trailing  edges  of  the  airfoil.  For  more  discussion 
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about  this  condition  refer  to  the  description  of  the  circular 
cylinder . 

Circular  Cylinder/ Compressible  Flow 

The  description  of  the  compressible  flow  problem  for  the  cylinder 
was  purposely  placed  after  that  of  the  airfoil.  The  problem  of  sym¬ 
metric,  steady,  potential  flow  over  a  circular  cylinder  can  be  formu¬ 
lated  in  half-space  (denoted  ,  as  described  for  the  incompressible 

problem.  The  governing  equation  when  compressibility  is  considered  is 
no  longer  the  Laplace  equation,  but  is  highly  nonlinear  and  of  mixed 
type  (Ref  24) .  The  exact  potential  equation  is 

+  2  y 

+  (#,y-al)i,,y  =  O 


where 

a'-.  &+  ^K-a:+i.y)] 

This  equation  is  a  nonlinear,  second -order ,  mixed  elliptic-parabolic- 
hyperbolic,  partial  differential  equation.  It  reduces  to  the  Laplace 
equation  for  incompressible  flow  and  to  the  small -disturbance  velocity 
potential  equation  (i.e.  eq  11)  for  flow  over  slender  bodies.  Instead 
of  using  the  exact  governing  equation  for  the  compressible  problem,  the 
small-disturbance  potential  equation  for  transonic  flow  was  selected  as 
the  governing  differential  equation.  From  a  physical  viewpoint  eq  11 
does  not  accurately  model  the  flow  over  the  cylinder  since  the  velocity 
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per turbat ions  become  large.  From  a  mathematical  perspective  it  can 
be  solved  for  the  cylinder  when  appropriate  boundary  conditions  are 
specified.  The  physical  validity  of  the  solution  can  be  evaluated  by 
comparison  with  solutions  of  the  exact  potential  equation.  Thus,  the 
assumed  governing  differential  equation  for  compressible  flow  over  the 
cylinder  is  given  by  eq  11 

+  (fcyy  =  O  <u> 

The  cylinder  radius  is  used  in  place  of  the  chord  length  to  nondimens ional- 
ize  eq  11  for  the  cylinder  problem. 

The  boundary  conditions  for  the  synmetric  compressible  problem  are 
described  in  a  previous  section.  They  are  given  by  eqs  3,  5,  and  10: 


^  (fr  — >  o 

as  r  — »  00 

(3) 

4>,r  +  CoS  & 

=  O  for  (f,  0)  in  iSXc. 

(5) 

<t>,6  -  ° 

for  (jC,  d)  in  ^-&c. 

(10) 
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Flowfield  Discretization 

A  finite  element  solution  of  the  governing  differential  equation 
for  the  cylinder  problem  requires  that  the  flowfield  be  discretized  by 
a  finite  number  of  elements.  Two  techniques  are  possible;  for  each 
technique  domain  SI.  is  replaced  by  a  finite  domain  Si  p  .  First  of  all , 
one  could  take  SLp  to  be  very  large,  and  require  that  the  actual 
gradient  boundary  conditions,  expressed  by  eq  3,  be  enforced  along  the 
farfield  boundary.  Since  these  boundary  conditions  and  also  those 
specified  on  the  airfoil  surface  are  of  the  Neumann  type,  then  the 
solution  of  the  governing  differential  equation  can  only  be  determined 
to  within  an  arbitrary  constant.  If  this  technique  was  used,  then  one 
must  specify  the  value  of  (j>  at  some  point,  preferably  at  a  nodal 
point  along  the  farfield  boundary,  so  the  solution  can  be  uniquely 
determined  everywhere. 

The  second  possible  technique  is  to  impose  along  the  farfield 
boundary  ^  hQ? )  the  condition  c  (fipp  The  expression  for  (ftp? 
should  be  an  asymptotic  solution  which  satisfies  the  infinity  con¬ 
dition  and  is  valid  in  the  farfield  of  SX  (Ref  18) .  This  approach 
is  commonly  adopted  by  investigators  using  finite  difference  methods 
(Refs  1 ,  10 ,  40) .  It  is  also  employed  by  others  using  finite  element 
methods  (Refs  9,  12,  13)  and  will  be  used  in  the  present  analysis. 
Generally  speaking,  a  relatively  smaller  domain  SL  p  is  required  for 
the  second  technique  than  for  the  first.  This  means  for  a  desired 
degree-of-accuracy  fewer  degrees -of-freedom  are  needed  to  solve  the 
problem,  which  translates  into  lower  computational  costs. 
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A  number  of  different  elements  can  be  used  to  discretize  domain 
St  p  .  When  triangles  (Ref  3)  or  quadrilaterals  (Ref  13)  are  used,  the 
shape  of  domain  Sip  is  altered  since  the  discretization  does  not 
exactly  describe  the  boundary  shape.  Ihis  approximation  may  not  introduce 
significant  error  in  the  solution  unless  large  elements  are  used  near 
the  boundaries.  The  isoparametric  quadrilaterals  described  by  Raju 
(Ref  27)  or  those  developed  in  this  study  will  exactly  discretize  domain 
SLf  .  For  this  reason,  the  element  shown  in  Fig  3  will  be  used  for 
the  cylinder  problems  considered.  Figures  4  and  5  show  the  discretiza¬ 
tion  of  Stf  for  problems  of  symmetric  flow  and  flow  with  circulation, 
respectively. 

Incompressible  Flow 

Before  describing  the  finite  element  analyses  of  this  problem,  a 
further  simplification  will  be  discussed. 

Superposition.  Mien  incompressible  flow  with  circulation  is  con¬ 
sidered,  the  velocity  potential  function  may  be  represented  by  two 
simpler  component  functions.  Since  the  governing  differential  equation 
and  boundary  conditions  are  linear,  then  the  perturbation  potential 
function  can  be  expressed  as 

4>  -  +  r<j>t  at) 


This  superposition  of  solutions  will  produce  two  problems,  each  of 
which  can  be  solved  without  concern  about  what  the  value  of  circulation 
should  be  or  Itow  it  should  be  handled  numerically.  The  function  4* 
is  the  solution  for  the  thickness  or  symmetric  problem,  and  the 
solution  of  the  lifting  problem  with  unit  circulation.  The  actual 


Figure  4  -  Discretization  of  Flowfield  for  Flow  Over 
a  Circular  Cylinder  without  Circulation 


value  of  circulation  tr)  is  determined  from  the  stagnation  condition, 
eq  7,  after  each  of  the  two  component  solutions  are  known. 

Substituting  eq  16  into  the  governing  differential  equation  and 
boundary  conditions  one  obtains  two  problems.  First,  the  s wine trie 
or  thickness  problem  is  defined  by 


DE: 

—  0  for  (] f!  &)  in 

(17) 

BC's: 

d.  f  toS  &  =  O  for  (\T|  0)  in 

(18) 

<kt6  =  0 

for  (r,  a)  in 

(19) 

=  4s 

for  (X,  G>)  in  ^Hfp 

(20) 

Secondly , 

the  lifting  problem  is  defined  by 

DE: 

vz4t(r> 8) 

—  O  for  in  -fLp 

(21) 

BC's: 

4>ji.r  -  o 

for  in  ^ Sic. 

(22) 

4n  ~ 

for  (r,  0)  in  dSlfp 

(23) 

K-  f 

-  1  for  (r  2  / ,  ©  =  0  * ) 

(24) 

The  functions  and  ^  are  farfield  expressions  for  the  thickness 
and  lifting  problems,  respectively.  They  represent  velocity  potentials 
for  a  source  and  vortex  of  unit  strength.  After  ^  and  A  are 
known,  P  is  determined  from  eq  7.  Substituting  eq  16  into  eq  7  and 


25 


solving  for  f1  gives 


r  - 


rs'ne  -  4>+te 


(25) 


'.e 


<s*)*0  ,-fc) 


Finite  Element  Solution.  Suppose  .Tip  is  discretized  by  a  total  of 
E  elements  with  a  total  of  |sj  system  nodes .  The  approximate  solution 
of  the  Laplace  equation,  which  governs  each  of  the  two  problems  described, 
is  obtained  from  the  method  of  weighted  residuals  as  expressed  by 

If  6)  ^  (fi  6)  cM  r  O  (26) 

Af 

for  i  :  l,  * ' N  •  Functions  4V  M)  are  weight  functions  that  will 
be  specified  later  by  the  Galerkin  method.  Integration  by  parts  or 
using  Green's  Theorem  gives 

Jj  dA  -  O  (27) 

SL  p 

The  functions  (fa  and  are  chosen  from  a  class  of  functions  so  that 
eq  27  is  integrable.  When  is  discretized  by  finite  elements, 

then  and  YV  must  be  continuous  across  inter-element  boundaries  and 
have  measurable  first  derivatives  throughout  St  p  ■  This  means  that 
(f)  and  cannot  be  piecewise  constants,  but  they  can  be  piecewise 
linear  or  bilinear  functions,  providing  continuity  is  enforced  across 
inter -element  boundaries. 

The  Galerkin  method  of  weighted  residuals  is  used  to  obtain  the 
finite  element  solution.  Piecewise  trial  functions  are  chosen  to 
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approximate  (f)  which  satisfy  both  the  continuity  requirements  and  the 
essential  boundary  conditions  (i.e.  boundary  conditions  for  which  is 
specified).  Within  each  sector  element,  denoted  £  ,  the  trial 
function  can  be  expressed  as 


4>%,e>)-  NjOr,©)^. 


(28) 


for  Js),*",  4  *  repeated  index  J  indicates  sunrnation  over  the 

range  of  the  index.  (\{j  ( ©)  are  basis  or  shape  functions  and  are 
the  unknown  nodal  values  of  the  potential  function.  For  Galerkin's 
method  the  weight  functions  ^  are  set  equal  to  the  shape  functions 
N;rte)  Since  the  basis  functions  are  chosen  to  be  continuous 
across  inter-element  boundaries,  then  eq  27  can  be  written  for  each 
element  as 


||  VIS/,  cM  -  j  ^4^NidS=0  (29) 

ne  an* 

for  t  ~  I ) • ‘  '  j  ^  The  star  -tf-  notation  on  the  boundary  term 
signifies  that  the  term  can  be  non-zero  only  for  elements  6.  that 
border  the  boundary  of  Sip  .  For  all  other  elements  the  boundary 
integral  is  zero.  Substituting  eq  28  into  eq  29,  and  also  putting 
the  actual  gradient  boundary  conditions  into  the  boundary  term,  pro¬ 
duces  elemental  equations  of  the  form 


(30) 


The  matrix 


is  referred  to  as  the  elemental  stiffness  matrix  and 
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is  given  by 

K;;  -  II (  N‘>  N;,r  +  7C  N;,s  HJ,s)jA  (31) 

Sle 

The  form  of  vector  -JT*  depends  on  the  gradient  boundary  conditions  for 
the  problem  being  considered  and  will  be  specified  in  later  sections. 

The  essential  boundary  conditions,  which  are  imposed  on  the  trial 
functions,  are  enforced  after  the  assembly  of  elemental  equations. 

Thus,  in  the  boundary  term  of  eq  27  is  selected  to  be  zero  along 
segments  of  c>XL  f  where  is  specified  (i.e.  the  farfield  boundary, 
ZSLff  ) .  The  usual  assembly  procedure  is  used  to  transform  all  of 
the  elemental  equations  into  a  global  system  of  equations  which  can 
be  expressed  as 

-  f;  ■■■'")  (32) 

Global  expressions  K  (fit  F*  are  the  counterparts  of  elemental 
expressions  X  f  q>  f  respectively.  Equation  32  is  reduced  to 
the  final  set  of  system  equations  by  enforcing  the  essential  boundary 
conditions  along  dXlpp  First,  eq  32  is  partitioned  as  follows 


- 1 

I 

:  Kat~ 

m 

M 

L»<u 

I 

!&. 

1  I 

IM 

(33) 


Vector  ^  (j>^  is  composed  of  the  L  nodal  values  of  which  lie 

along  &£2.pf  and  is  computed  from  the  farfield  expression 

ff 

Vector  ^  is  composed  of  the  remaining  N~  L  unknown  nodal 
values  of  (ft  ,  which  is  determined  from  eq  33  by  inverting  matrix 
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[K„]  t°  8^ 

[4]  =  [k«»]  (m_  a*) 

The  assembly  and  reduction  procedure  just  described  will  be  followed 
throughout  this  report. 

Trial  Functions.  Three  trial  functions  were  chosen  to  approximate 
the  solution  of  within  the  sector  element  shown  in  Fig  3-  They 


;e 

(i)  cf>C r}e)~  0,  -f  |>, 


~r-  +  o,  +  d,  e 


(2)  =  Q2+  ba.r  +  0 -b  dz re  (36 

(3)  0C(v-,©)  r  (X3  +  b3r'  +  03e  +  d3  re  (37 

Constants  (  Oi^  '**  dJ*  )  can  expressed  as  functions  of  the  unknown 
nodal  value  of  (jf>  and  the  geometric  parameters  of  the  element. 

When  the  constants  are  evaluated,  each  of  the  trial  functions  can  be 
written  in  the  form  expressed  by  eq  28 

Nj<n.)£/  (j=  <28 

The  expressions  for  N*  (^  ©)  are  given  in  Appendices  A-C  for  each 

j 

of  the  above  trial  functions.  Throughout  the  remainder  of  this  report 
the  superscript  G.  will  be  omitted  when  it  is  clearly  understood  that 
elemental  quantities  are  being  considered. 

Symnetric  Flow.  First,  consider  the  problem  of  synmetric  flow. 
This  problem  can  be  formulated  by  either  considering  all  of  domain 
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rzf  or  half  of  it.  Typical  discretizations  are  shown  in  Figs  U 
and  5.  Using  Galerkin's  method  eq  27,  as  expressed  by  eq  29.  can  be 
written  for  element  £.  as 


4 


II  (  ^  N;  •  7  Nj  )  a  A  4> >  •  =  I  V<l>>n  N;  dS  (38) 

for  The  star  notation  on  the  boundary  term  sig¬ 

nifies  that  the  term  is  present  only  for  elements  which  intersect 

and  When  eq  18,  the  boundary  condition,  is  substituted 

into  eq  38,  the  results  can  be  written  as 


(30) 


where  stiffness  matrix  is  given  by  eq  31.  The  only  non-zero 

ne  J 

values  for  +,'  canes  from  elements  which  contact  the  cylinder  surface 
and  is  defined  by 

^ N;  Cos>&  rjc)0  ( 


The  symnetric  problem  was  solved  using  the  three  sector  elements 

e  re 

described  in  Apendices  A-C.  The  derivations  of  and  t.  are  also 

4 

presented  in  these  appendices.  The  assembly,  reduction,  and  solution 
procedure  described  in  the  previous  section  was  followed.  The 
expression  for  the  farfield  potential,  needed  to  reduce  the  assembled 
equations,  is  given  by 
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<t>Ff  --  4  --  tosA 


for  (X,BS)  in 


(40) 


where  v*  Rff  is  the  radius  of 

Flow  With  Circulation.  When  circulation  or  lift  is  present  the 
entire  domain  Slf  must  be  discretized  as  shown  in  Fig  5.  Along  the 
branch  cut,  nodes  are  placed  at  pair-points  located  on  opposite  sides 
of  the  cut.  The  difference  in  potential  between  pair-point  nodes  must 
be  equal  to  the  value  of  circulation. 

A  convenient  way  to  consider  flow  with  circulation  is  to  use  the 
superposition  technique  previously  described.  The  thickness  or  sym¬ 
metric  problem  has  already  been  discussed.  The  lifting  problem, 

described  by  eqs  20-24,  is  formulated  numerically  in  a  similar  manner. 

e 

The  elemental  stiffness  matrix  is  the  same  as  for  the  thickness 

problem,  since  the  governing  differential  equation  for  both  problems 
is  the  Laplace  equation. 

There  are  four  minor  formulation  differences  between  the  two 
problems.  First,  the  vector  equals  zero  for  every  element  £ 
for  the  lifting  problem.  This  is  due  to  the  boundary  condition  speci¬ 
fied  along  ,  eq  22.  Secondly,  the  farfield  expression  is  given 

by 

for  (T)  ^  )  in  ZSlff  (41) 

Thirdly,  the  required  lump  in  potential  across  the  branch  cut  must  be 
enforced.  This  is  done  bv  setting  the  nodal  values  of  z.  O  and 
c  -  I  for  nodes  along  the  cut  This  choice  is  consistent  with 
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the  farfield  expression  as  y0  — >  &FF  and  satisfies  the  requirement 
specified  by  eq  24.  This  choice  also  forces  c  O  for  points 

along  the  cut.  Since  this  condition  is  true  everywhere  for  purely 
circulatory  flow,  then  no  generality  is  lost  by  its  use.  It  should 
be  pointed  out,  that  this  choice  would  not  be  possible  if  the  super¬ 
position  principle  were  not  used  to  divide  the  total  problem  into  two 
simple  problems.  When  the  total  problem  is  formulated  without  super¬ 
position,  then  nodal  values  of  potential  along  both  the  upper  or  lower 
sides  of  the  cut  must  be  left  unspecified.  This  will  insure  that  the 

t  r*  is  not  forced  to  be  zero,  which  it  is  not  for  the  thickness 
problem.  The  fourth  difference  occurs  in  the  reduction  procedure. 

Since  the  nodal  values  are  specified  along  the  cut,  then  they  have  to  be 
included  with  the  nodes  along  the  farfield  boundary  when  the  assembled 
equations  are  reduced.  Once  solutions  for  ^  and  ^  are  determined, 
then  the  circulation  defined  by  eq  25  can  be  calculated. 

The  problem  of  flow  with  circulation  was  solved  using  sector 
elements  (2)  and  (3) .  Appendices  B  and  C  present  derivations  of  the 
elemental  equations. 

Compressible  Symmetric  Flow 

The  problem  of  compressible  flow  without  circulation  is  formu¬ 
lated  in  the  half-space  shown  in  Fig  4.  The  governing  differential 
equation  for  this  problem  is  given  by  eq  11,  and  the  boundary  con¬ 
ditions  by  eqs  3,  5,  and  10. 

Discretization.  The  solution  of  this  problem  by  the  finite 
element  method  is  obtained  in  the  finite  half-space  XI  p  ,  as  described 
for  the  incompressible  problem.  Domain  is  discretized  as  shown 
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in  Fig  4  by  a  total  of  £  elements  with  fsj  global  nodes.  Since  the 
governing  equation  for  this  problem  is  the  same  as  the  equation  for 
the  airfoil  problem,  then  the  asymptotic  solution  used  as  a  boundary 
condition  along  the  far  field  boundary  (SSLpp  )  is  given  by  eq  50 
described  in  Chapter  IV. 

Finite  Element  Solution.  From  the  method  of  weighted  residuals, 
eq  11  can  be  expressed  as 

//{[ 1  ',vV i*]  +  <t>»]  1>-  d  »  =  o  («) 

Slf 

for  cz  ),  • " ,N-  For  appropriately  chosen  weight  functions  this  ex¬ 
pression  can  be  integrated  by  parts  to  give 

IK  -  <t>X )K,*]  J a 

©  (M) 

ail*  J 

The  finite  element  approximation  for  in  each  element  is  given  by 
eq  28  as 

4>(y,  6)  -  Nj  (r,e>)  (f>j  (28) 

for  j  -  For  Galerkin’s  method  the  weight  functions  arG 

set  equal  to  the  shape  functions  fsjj .  When  shape  functions  are 
chosen  to  insure  integrability  of  eq  43,  then  it  can  be  written  for 
every  element  and  assembled  as  previously  described.  The  elemental  form 
of  eq  43  can  be  written  in  iterative  form  and  expressed  symbolically  as 


(44) 


[ttj  +  L;j  (j*)J  4.  =  f’  +  Jt'(0 


for  ((j's  I,  4  The  elemental  matrix  and  vector  p;  are  the 

same  as  those  defined  for  the  incompressible  synmetric  problem  The 

new  expressions  and  come  from  the  nonlinear  term 

J  ** 

and  are  functions  of  the  solution  itself.  The  superscript  V\  is 
associated  with  the  iterative  solution  procedure  which  is  briefly 
described  below  and  discussed  in  more  detail  in  Chapter  IV.  The  super - 
script  Y\  should  not  be  confused  with  the  vector  y\  which  is  the  unit 
outward  normal  along  the  cylinder  surface. 

The  iterative  form  of  the  elemental  equations  comes  from  the  way 

the  governing  equation  is  solved.  The  nonlinear  term,  *  in  eq  43, 

,v\  ,t\  +  l 

is  written  iteratively  as  0*  ©  Essentially,  the  potential 

A  “/  XX 

function  is  replaced  by  a  sequence  of  functions 

r.r'i  which  generate  a  sequence  of  equations  expressed  by  eq  44. 
These  equations  are  solved  for  each  iteration  until  the  sequence  of 
potential  functions  converges.  Convergence  is  assumed  when 


<t>* 


f:  e 


(45) 


for  some  small  €.  . 

An  additional  approximation  is  made  to  simplify  the  integrations 
required  by  eq  43.  In  the  nonlinear  term  the  derivative  Is  as- 

sumed  to  be  a  constant  in  each  element.  The  constant  chosen  is  the 
average  value  of  <pt%  computed  from 
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This  approximation  is  made  only  for  the  cylinder  problem.  The  non¬ 
linear  term  is  not  locally  linearized  in  this  sense  for  the  airfoil 


problem,  but  is  treated  in  a  more  exact  manner.  For  this  approximation 


V 

matrix  (.  l  •  (<*")  can  be  expressed  in  two  parts  as  +" 

c  ■  J  '  J  ^ 

^;C4>y  These  expressions  are  given  by 


L;.  (*-)  c  -  n4,  [|+  k2  0  +  *)J  JJ  [C“  6  N‘,r  W,',r 

J  n  ~ 


lie 


5ifl  (9 

r 


& 


ra“(N*>'&»  +  *>»! i>) 

N‘>Nj',eJ 


(47) 


and 


LyW*)  =  -  N- n.i6  r| 

Vector  is  given  by 


r=  t 


(48) 

(49) 


rsijsl 


^  i  - 

The  synmetric  compressible  proble.  was  solved  using  sector  element 


(3).  Appendix  C  contains  a  description  of  the  element  and  the  deriva¬ 
tives  of  the  elemental  equations. 
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IV  Analysis  of  Flow  Over  a  Thin  Airfoil 


The  governing  equation  and  boundary  conditions  for  flow  over  a 
thin  airfoil  are  described  in  Chapter  II.  This  chapter  describes  the 
numerical  analysis  of  the  problem  by  the  Finite  Element  Method.  Many 
of  the  formulation  techniques  and  procedures  that  are  used,  have  been 
described  in  Chapter  III  for  the  cylinder  problem.  Similar  discussions 
are  not  repeated  for  the  airfoil;  instead,  the  reader  should  refer  to 
the  appropriate  sections  of  that  chapter  for  more  detail. 


Flowfield  Discretization 

The  infinite  flow  domain  n.  is  replaced  by  a  finite  but  suffi¬ 
ciently  large  dcmain  I"Lp  with  the  potential  specified  along  the  far- 
field  boundary  d  pe  .  A  detailed  discussion  of  the  reasons  for  this 
approximation  is  given  in  Chapter  III.  The  condition  imposed  along 
^  -Tb  f  p  is  the  farfield  expression  of  Klunker  (Ref  18)  given  by 


fF  "  TTf,  + 


r[ 


chord 


pS~s}"(9) 


+  j'An'  ~~~ 


J±L(f  - j,  u 


where  \ 
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The  first  term  in  eq  50  is  directly  associated  with  the  thickness 
distribution  (  *C  is  the  thickness  ratio)  of  the  airfoil,  and  is  the 
dominant  term  for  nonlifting  flow  (Ref  18) .  For  lifting  flow,  the 
second  term  is  the  most  dominant  of  the  three  and  was  the  only  term 
included  in  Chan's  formulation  (Ref  13).  The  expression  sgn(y')  is 
required  to  account  for  the  jump  of  potential  in  the  farfield  down¬ 
stream  of  the  trailing  edge.  For  lifting  flow,  a  branch  cut  must  be 
included  in  Si  p  extending  from  the  trailing  edge  to  the  farfield 
boundary.  Across  this  cut  the  jump  of  potential  is  forced  to  equal 
the  value  of  F  appearing  in  eq  50.  The  last  term  is  of  highest  order 
and  involves  an  area  integration  over  the  entire  domain  SLf  .  When 
this  term  is  included  in  the  farfield  expression,  an  iterative  solution 
algorithm  must  be  used  since  the  term  depends  upon  the  solution  itself. 
The  effect  of  neglecting  the  last  term  is  evaluated  in  Chapter  VI. 

Consider  the  problem  of  flow  over  a  symmetric  airfoil  at  zero  angle 
of  attack.  From  symmetry  considerations  the  problem  can  be  formulated 
in  the  half-space.  Figure  6  shows  an  initial  division  of  into 

three  segments  or  super-elements.  Since  the  tangential  boundary  con¬ 
dition  for  a  thin  airfoil  is  enforced  along  the  ©  axis,  then  the 
first  row  of  elements  used  to  discretize  the  center  segment  must 
extend  to  the  chord  line  (  %-  axis).  Any  type  of  element  with  straight 
line  boundaries  (i.e.  triangles,  rectangles,  etc.)  could  be  used  to 
discretize  the  super-elements.  After  selecting  the  type  of  element 
to  be  used,  one  could  write  a  computer  subroutine  to  automatically 
discretize  each  super-element  from  a  few  input  parameters  which  describe 
how  the  segment  is  to  be  divided  into  smaller  parts.  This  procedure  was 
followed  for  both  the  cylinder  and  the  airfoil  problem;  however,  the 
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details  are  not  important  and  will  not  be  discussed. 

A  slight  modification  in  the  described  discretization  of  the 
center  super-element  would  be  required  if  the  tangential  boundary  con¬ 
dition  was  applied  along  the  actual  airfoil  contour.  For  this  case, 
elements  in  the  center  segment  would  not  extend  to  ^cO  .  Elements 
would  have  to  be  selected  to  either  match  or  approximately  match  the 
shape  of  the  thickness  profile.  This  difference  in  discretizations 
becomes  significant  when  thick  airfoils  are  to  be  analyzed,  and  is 
discussed  further  after  the  elemental  equations  are  formulated. 


Iterative  Approximation 

If  the  solution  of  eq  11,  the  governing  differential  equation,  is 
directly  formulated  by  finite  element  methods,  then  a  set  of  second- 
order,  nonlinear,  algebraic  equations  will  result.  These  equations 
would  have  to  be  solved  by  some  iterative  technique  such  as  Newton- 
Raphson.  This  type  of  solution  process  can  be  avoided  by  directly 
expressing  and  solving  the  governing  equation  in  iterative  form.  First, 
eq  11  is  written  as 

«)<*,*  -  it!  +  i jtyy  =  O  (51) 

As  described  previously  for  the  cylinder,  the  nonlinear  term  ^  is  ap- 

jY\  +  t 

proximated  iteratively  by  (pf  ^  (p*  ^  where  f\  denotes  the  iteration 
number.  In  essence  potential  ^  is  replaced  by  a  sequence  of 
potentials  \  4>°,  4 \  ,  which  converges  when 
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for  small  £  .  Thus,  eq  51  written  in  iterative  form  becomes 


r  T~  v  +  >  r  iiv  /*1  in  +  ' 

ii  J  +  ^yy  =  o  (52) 


Finite  Element  Solution 

Suppose  the  half-space  SLp  is  discretized  by  a  total  of  £ 
elements  with  N  nodes.  The  approximate  solution  of  eq  52  is  ob¬ 
tained  from  the  method  of  weighted  residuals  as  expressed  by 


SLf 


l±r 

2. 


A  =  o 

(53) 


where  4  x  I;'*’, N  •  For  suitable  weight  functions  'fc1  .  this  equation 
can  be  integrated  by  parts  to  give 


ua  ya  I  ^  ^  I  w 


a.a« 


(54) 


where  y\  -  (  Y)  is  the  unit  normal  vector  along  each  segment  of 

The  boundary  conditions: 


<t>ty  -  (l  +  )'3ym  for  in  (55> 


r»-H 
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and 


for  (X<V)  in  ^-^-o 


(56) 


are  substituted  into  the  boundary  term  of  eq  54.  Since  <f)  is  specified 
along  ,  then  w  is  taken  to  be  zero  there. 

Equation  54  is  piecewise  integrable  over  the  discretized  half¬ 
space  providing  and  are  at  least  continuous  functions  across  the 
inter- element  boundaries.  Within  each  element  the  solution  can  be 
approximated  by 


4>(*- y)  -  Nj  (x,v)  <(>j 


(57) 


where  |S|  •  are  the  shape  functions  and  q>  •  the  unknown  potential  values 
at  the  nodes.  The  shape  functions  are  chosen  to  satisfy  the  required 
continuity  of  ^  in  the  global  sense.  For  Galerkin’s  method  is  set 
equal  to  Thus,  all  continuity  requirements  are  satisfied  to 

allow  eq  54  to  be  written  in  elemental  form  and  assembled  to  obtain 
the  global  form.  The  elemental  form  of  eq  54  is  expressed  as 

K(j  (<fn)  4>J*'  -  t  <58> 

LX  Km  is  given  by  the  sum  of  the  following  five 

J 


Matrix 


matrices : 


Aa  -  (  h  [j*  dxJy  (58a) 

J  SU 


Mj;y  d  X  dy  (58b) 
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^<j  ($*)  -  II  cl/dy 

J  Sip 


(58c) 


fcj  =  <  bN;  Vl£|  dx 

A  \l  -  /-i 


(58d) 


\|=o 


Etj(f)  =  Mi  ^  £  I  N*<*  N/  N',^1  Jx 

5< 


(58e) 


Vector  *pj  is  defined  by 

t  =  !,  N.'#l  i* 


(58f) 


iXvA  y=  o’ 


In  general  |^«  (**)  and  •£'  depend  upon  the  shape  of  the  airfoil. 

J  ' 

The  contribution  to  k,v  m  from  the  airfoil  shape  comes  only  from 

J 

matrices  DC‘J  and  ^M")  ■  These  matrices  are  evaluated  only 
for  the  first  row  of  elements  above  the  airfoil  in  the  center  segment 
shown  in  Fig  6.  They  are  zero  for  all  other  elements  since  they  come 
from  the  global  boundary  term  of  eq  54.  This  boundary  term  could 
alternatively  be  taken  to  the  right-hand  side  of  the  equation  and 
treated  as  a  force.  The  effect  of  this  alteration  would  be  a  slightly 
slower  convergence  rate  caused  by  replacing  +  >  with  <$> *  in  the 
affected  terms. 

The  contribution  to  KtV^*)  from  matrices  A.'; ,  8.’;  and 

j  /  j  '  j  * 

are  independent  of  the  airfoil  being  considered,  since  they 

A 

depend  on  an  integration  over  the  elemental  area.  This  is  true  only 
because  the  tangential  boundary  condition  is  imposed  along  ^  c  O  ^  , 
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which  requires  that  elements  be  extended  to  the  X  -axis.  If  the 
boundary  condition  was  enforced  on  the  airfoil  surface,  then  the  area 
of  elements  in  the  first  row  of  the  center  segment  would  be  smaller 
by  comparison.  Thus,  an  indirect  dependence  of  the  airfoil  shape  from 
these  elements  would  be  reflected  by  all  the  matrices  defining 
K;;0n  The  effect  of  enforcing  the  boundary  condition  along 
^  -  O*  instead  of  on  the  airfoil  contour,  is  not  critical  for  thin 
airfoils  as  long  as  element  sizes  in  the  ^-direction  are  larger  than 
the  thickness  of  the  airfoil.  Since  one  of  the  advantages  of  the  FEM 
is  the  ability  to  use  relatively  large  elements  to  achieve  accurate 
approximations,  then  element  size  difficulties  should  not  be  critical 
until  thickness  ratios  become  large.  For  thick  airfoils  or  for  arbi¬ 
trarily  shaped  bodies  the  boundary  conditions  must  be  satisfied,  not 
on  ^  -  O  ^  >  but  on  the  surface  contour .  For  these  shapes ,  elements 
would  terminate  at  the  contour  and  pose  no  conflict.  By  way  of  com¬ 
parison,  some  investigators  (Refs  10,  40)  using  finite  difference 
methods  and  small  grid  sizes  satisfy  the  tangential  condition  along 

'j-  ot 


Several  different  types  of  elements  and  orders  of  approximation 
could  be  used  to  solve  this  problem.  For  example,  the  higher-order 
cubic  triangular  elements  used  by  Chan  (Ref  13)  would  be  more  than 
sufficient.  Also,  linear  triangular  elements  (Ref  7),  which  provide 
the  lowest  permissible  approximation,  would  be  adequate  since  they 
satisfy  the  required  continuity.  In  terms  of  velocities,  the  linear 
triangles  would  give  constant  velocity  elements.  In  order  to  get  a 
linear  variation  in  velocity,  which  is  not  required,  one  would  have  to 
use  quadratic  triangular  elements.  Somewhere  between  the  linear  and 
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quadratic  triangular  elements  are  the  bilinear  rectangles  or,  in 
general,  the  quadrilaterals.  Rectangles  are  normally  not  very  useful 
except  for  simple  boundary  geometries. 

For  the  problem  being  considered,  the  boundaries  are  straight  line 
segments  as  shown  in  Fig  6.  The  boundary  segment  between  the  leading 
and  trailing  edges  of  the  airfoil  is  straight  only  because  the  tangen¬ 
tial  boundary  condition  is  satisfied  along  the  chord  (i.e.  y  -z  o  +  )  • 
Therefore,  rectangles  can  be  used  effectively  to  discretize  -Qf 
everywhere.  Rectangles  are  desirable  due  to  their  simplicity,  and  were 
chosen  for  this  reason.  They  are  described  in  Appendix  D  along  with 
the  derivation  of  elemental  equations. 

Mixed  Flow 

Reported  Convergence  Difficulties.  Several  investigators,  who 
have  tried  to  solve  the  velocity  potential  equation  using  finite  element 
methods,  report  difficulties  with  convergence  of  solution  algorithms 
for  transonic  flow.  Among  these  are:  Chen  and  Habashi  (Ref  12),  Chan 
(Refs  13-16) ,  Ecer  (Ref  22) ,  Akay  (Ref  24) ,  and  Aziz  (Ref  23) .  Diver¬ 
gence  of  solution  algorithms  have  lead  some  to  believe  that  Galerkinfs 
method  could  not  be  used.  Others  claim  the  source  of  difficulty  is 
the  small-disturbance,  velocity  potential  formulation  of  the  problem 
and  suggest  that  alternative  formulations  be  tried.  However,  Akay 
(Ref  24)  reports  convergence  difficulties  with  finite  element  solutions 
of  the  total  velocity  potential  equation.  It  should  be  noted  that  con¬ 
vergence  problems  also  occur  when  finite  difference  methods  are  used 
(Refs  33,  34,  40).  To  insure  that  solution  algorithms  converge, 
special  difference  operators  have  been  developed.  A  different  operator 
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may  be  used  at  each  grid  point  depending  on  whether  the  point  is  con¬ 
sidered  to  be  an  elliptic,  parabolic,  hyperbolic,  or  shock  point 
Divergence  of  solution  techniques  for  more  classical  methods,  such 
as  the  Rayleigh- Janzen  Method  (Refs  29 ,  30) ,  also  occurs  when  transonic 
flow  develops. 

The  convergence  problems  reported  by  other  investigators  for 
transonic  flow  are  also  observed  when  the  methods  described  in  this 
study  are  used  to  solve  the  transonic  problem.  Solution  techniques 
converge  quickly  as  long  as  the  flow  remains  subsonic  everywhere  in 
Sic  ■  When  transonic  flow  occurs  (i.e.  a  small  supersonic  bubble 
appears  in  the  flow) ,  solution  techniques  which  are  suitable  for  sub¬ 
sonic  flow  (i.e.  an  elliptic  problem)  do  not  converge  at  all.  These 
conclusions  are  supported  by  solution  results  presented  and  discussed 
in  Chapter  VI.  Convergence  difficulties  are  caused  by  the  mixed 
character  of  the  flow  which  results  when  the  governing  differential 
equation  changes  type  from  elliptic  to  hyperbolic  in  the  supersonic 
region.  To  overcome  convergence  difficulties  an  adjustment  is  needed 
in  the  finite  element  equations  for  elements  in  the  supersonic  region. 
The  purpose  for  this  adjustment  is  to  account  for  the  proper  zones  of 
influence  in  the  supersonic  region  where  the  equation  is  hyperbolic. 
Finite  element  formulations  of  elliptic  equations  work  quite  well,  but 
those  same  formulations  applied  to  hyperbolic  equations  will  not  work. 

Formulation  Adjustments,  Other  investigators  have  tried  various 
techniques  to  alter  their  finite  element  formulations  for  mixed  or 
transonic  flow.  These  techniques  are  briefly  described,  and  all  were 
tried  in  the  present  study.  The  first  technique  was  developed  by 
Chan,  Brashears,  and  Young  (Ref  23).  They  altered  their  finite  element 
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equations  for  elements  in  the  supersonic  zone  during  the  assembly 
procedure.  Before  the  local  stiffness  matrix  was  assembled  the 

coefficient ,  (J  £  (H'TTMk  ,  was  calculated 

at  each  node.  If  C,  was  negative  at  all  nodes  in  the  element,  then 
the  rows  in  the  stiffness  matrix  corresponding  to  the  downwind  in¬ 
fluence  were  ignored  in  the  assembly  process.  By  zeroing  out  the 
appropriate  rows  in  the  elemental  stiffness  matrix  the  downwind  in¬ 
fluence  on  the  solution  at  upwind  nodes  was  blocked.  Chan’s  technique 
appeared  to  work  in  conjunction  with  the  least  squares  formulation 
that  he  used. 

The  second  technique  was  developed  by  Akay,  Ecer,  and  Utku  (Ref 
24).  For  the  iterative  solution  procedure  they  chose,  the  stiffness 
matrix  was  a  function  of  the  nerturbation  velocity  squared.  During 

the  iterative  process,  if  any  alemer  Q  was  inside  the  supersonic 

Z 

zone,  then  the  velocity  term  denoted  by  was  replaced  by 

(I"  (59) 

Velocity  is  the  velocity  in  the  element  on  the  upstream  side  of 
element  £  .  The  constant  9  is  the  upwinding  coefficient  which  was 

taken  to  be  between  0.2  and  0.3.  According  to  Akav,  this  technique 
prevented  the  iterative  solution  algorithm  from  diverging;  however,  it 
did  not  produce  convergence  in  the  sense  of  eq  45.  Instead,  the 
solution  for  the  potential  function  oscillated  about  some  solution. 

In  addition  to  Akay's  own  technique,  he  also  tried  Chan's  upwinding 
method,  which  caused  the  solution  algorithm  to  diverge  inmediately 
upon  application. 


r  rfr  ■  iT 


Four  additional  upwinding  techniques  were  tried,  three  of  which 
are  briefly  described  below.  The  first  is  an  element  sliding  tech¬ 
nique  which  is  implemented  during  the  assembly  process.  As  the 
elements  are  assembled,  a  check  is  made  to  see  if  the  velocity  at  all 
element  nodes  is  supersonic.  If  the  element  is  inside  the  supersonic 
region,  then  it  is  slid  upstream  before  assembly.  Thus,  the  influence 
on  the  unknown  nodal  values  associated  with  that  element  comes  frcm  the 
region  upwind  of  the  nodes. 

The  second  technique  is  a  nodal  sliding  method  which  is  imple¬ 
mented  after  partial  assembly  of  the  elements  in  the  supersonic  bubble. 
This  method  closely  resembles  the  upwinding  or  backward-difference 
methods  used  in  finite  difference  analyses.  Before  assembling  any 
elements,  a  check  is  made  to  identify  which  nodes  are  contained  within 
the  supersonic  bubble.  If  a  node  is  in  the  bubble  then  all  elements 
which  have  that  node  in  cannon  are  identified.  Next,  a  portion  of  the 
stiffness  matrix,  which  comes  from  the  X  -derivative  terms,  is  assem¬ 
bled  for  the  identified  elements.  After  this  portion  of  the  assembly 
process  is  conpleted,  then  the  coefficients  in  the  partially  assembled 
global  equations  are  slid  upwind  one  position.  This  procedure  is  anal¬ 
ogous  to  taking  the  first-order,  five-point,  central -difference 
star  and  sliding  the  horizontal  coefficients  upwind  one  step  to  get 
the  backward-difference  operator  as  illustrated  below 

*  # 

O  *  •  •  #  »  #  o 

•  * 
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The  third  upwinding  technique  was  only  briefly  examined.  It 
was  taken  from  a  suggestion  by  Christie  (Ref  37) ,  and  could  prove  to 
be  effective  if  developed  properly.  For  this  technique,  the  weight 
functions  are  chosen  to  be  different  than  the  shape  functions  for  all 
elements  inside  the  supersonic  bubble.  The  idea  is  to  give  more  weight 
(i.e.  more  influence)  to  the  upwind  half  of  the  element  than  to  the 
downwind  half.  For  the  technique  tried,  upwinding  was  done  only  in  the 
jC  -direction.  Weight  functions  were  taken  to  be  4*/  c.  |\|t*  +  o(  'JD(X) 
for  c  equal  to  an  upwind  node,  and  4V  ~  *“  for 

f 

C  equal  to  a  downwind  node.  The  function  PC*)  is  chosen  to  be  zero 
at  the  nodes.  This  technique  was  tried  for  7>(*)  =  x(x-  2  d) 
and  c4  x  I  ,  but  it  was  not  explored  thoroughly  and  warrants  further 
study. 

New  Upwinding  Technique.  A  new  upwinding  technique  was  developed 
in  this  study.  It  is  more  intuitive  than  analytical  in  nature,  although 
it  has  an  analytical  foundation.  It  is  not  as  elegant  an  idea  as  the 
previous  one,  but  it  is  simple  to  use,  and  provides  accurate  approxi¬ 
mations  of  pressure  distributions  for  transonic  flows.  For  this  and 
other  reasons  it  may  be  preferred  to  the  previous  method.  It  is  well 
known  that  picking  weight  functions  which  are  different  from  the  shape 
functions  may  lead  to  significant  error  (Ref  37)  . 

The  new  method  was  developed  by  modifying  the  finite  element  formu¬ 
lation  of  the  nonlinear  term,  which  is  the  term  that  is  responsible  for 
the  mi; r*  character  of  the  flow  for  subsonic,  free-stream  Mach  numbers. 

Recall  that  the  nonlinear  term  in  eq  51  was  written  iteratively 
as  €  €'  For  iteration  (n  +  i)  this  term  is  relaxed  by 
replacing  it  with  the  expression 
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(60) 


The  relaxation  coefficient  R  takes  the  range  of  values :  O  —  R  —  I  . 
The  last  term  in  eq  60  can  be  taken  to  the  right-hand  side  of  the 
governing  equation  and  treated  as  a  force  computed  from  the  previous 
iteration.  Thus,  eq  51  written  in  iterative  form  becomes 

When  R -  |  ,  the  nonlinear  term  is  not  relaxed  at  all.  For  Rr  O  , 
the  nonlinear  term  is  totally  relaxed  (i.e.  it  is  taken  to  the  right- 
hand  side  of  the  equation  and  treated  as  a  force  that  depends  upon  the 
previous  iteration  of  the  solution).  For  values  of  R  between  0>  O  and 
l,  0  a  mixture  cf  these  two  extremes  exists.  The  sole  modification 
of  relaxing  the  nonlinear  term  is  not  sufficient  to  make  the  solution 
algorithm  converge  in  the  sense  of  eq  45.  It  would  be  difficult  to 
imagine  that  a  hyperbolic  equation  could  be  solved  by  merely  solving  a 
sequence  of  elliptic  (Poisson)  equations.  The  ideas  of  domain  of 
dependence  and  range  of  influence  from  the  theory  of  differential 
equations  must  be  included  into  the  modification  process.  As  dis¬ 
cussed  previously  the  sign  of  the  A XX  coefficient  will  determine 
whether  the  governing  equation  is  elliptic  or  hyperbolic.  Inside  the 
supersonic  bubble  that  appears  for  transonic  flow  the  coefficient  is 
negative  and  the  equation  is  hyperbolic.  Outside  of  this  bubble  the 
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coefficient  is  positive  and  the  equation  is  elliptic.  For  elliptic 
differential  equations  the  domain  of  dependence  is  the  entire  domain 
of  the  problem,  but  for  hyperbolic  equations  the  domain  of  dependence 
at  some  given  point  is  confined  to  the  region  within  the  backward 
cone  defined  by  the  characteristic  curves.  Consequently,  the  solution 
at  the  given  point  is  influenced  only  by  the  solution  at  points  in 
the  backward  cone  and  not  by  the  solution  at  any  other  points  in  the 
domain.  Likewise,  the  solution  at  the  given  point  will  influence  the 
solution  only  for  points  in  the  forward  cone,  also  defined  by  the 
characteristic  curves. 

Without  getting  involved  with  the  method  of  characteristics,  the 
above  described  concepts  can  be  incorporated  into  the  modification 
procedure  being  discussed.  For  elements  inside  the  supersonic  bubble 
the  nonlinear  term  is  relaxed  as  described.  This  means  that  forces 
idoich  depend  upon  the  solution  are  applied  at  each  of  the  element 
nodes.  However,  the  solutions  at  upwind  nodes  cannot  be  influenced 
by  the  solution  at  downwind  nodes;  therefore,  the  forces  acting  at 
downwind  nodes  are  set  equal  to  zero.  In  addition,  the  solution  in 
the  element  cannot  depend  on  forces  applied  at  upwind  nodes,  when 
those  forces  are  determined  from  an  integration  over  the  entire  area  of 
the  element.  At  any  point  in  the  element  the  solution  should  depend 
on  the  potential  only  in  a  backward  cone  defined  by  the  characteristic 
curves.  Thus,  the  upwind  forces  must  be  reduced  in  magnitude  by  some 
factor.  One  way  to  represent  this  numerically  is  to  multiply  the 
relaxed  term  R  )  W* )  by  the  factor  U  where  0.0  £r  U  —  1.0 
at  upwind  nodes.  Thus,  the  iterative  Galerkin  form  of  the  governing 
equation  for  hyperbolic  elements  becomes 
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where  U  is  given  by 


U  ■ 


Q,0  to  | ,  O  at  upwind  nodes 
Q,  O  at  downwind  nodes 


In  actual  implementation  eq  62  is  integrated  by  parts  as  discussed 
previously.  The  elemental  equations  expressed  by  eqs  58  a-f,  and 
referred  to  as  equations  for  elliptic  elements,  are  modified  by  the 
upwinding  procedure  for  hyperbolic  elements  as  follows: 

(AcjKyp  -  (A‘j  )etL 

(8ij)Wyp  -  (&\j)eu 

Cy  tfJKyP=  R^)«u 

(Dij  )  kyp  ~  )eu 
Eij  Wkyp* 

CWkvp  =  tf‘)cu  t  (V+KOkvp 


( 


The  last  two  terms  in  eq  63f  are  forces  resulting  from  the  upwinding 
procedure.  They  are  given  by 


(?<)hyp  “ 

-  u Cl-  r)  <t>- 

(63p' 

-  0(|- R)  E.-jWeu 

(63h) 

The  relaxation  and  upwinding  parameters  IS.u)  were  investigated 
by  numerical  experiment  for  two  reasons.  First,  an  estimate  was 
needed  on  the  range  of  values  for  which  solution  algorithms  converge. 
Secondly,  the  "best  possible"  combination  of  parameter  values  for  a 
given  Mach  number  was  needed  to  solve  for  realistic  pressure  distri¬ 
bution  on  the  airfoil.  The  effect  these  parameters  have  on  conver¬ 
gence  properties  and  on  airfoil  pressure  distributions  is  presented 
in  Chapter  VI. 

In  general,  the  idea  behind  the  upwinding  method  is  to  modify 
the  formulation  just  enough  to  capture  the  physics  of  the  problem, 
within  the  realm  of  the  assumptions. 

For  a  given  airfoil  shape  and  Mach  number,  the  values  of  R  and 
yj  are  selected  by  an  iterative  method.  First,  upper  bounds  (i.e. 
values  above  which  solutions  will  diverge)  and  lower  bounds  (i.e. 
values  below  which  solutions  will  not  diverge)  are  estimated.  The 
idea  is  then  to  determine  the  lowest  upper  bound  for  both  R  and 
U  which  permits  the  solution  algorithm  to  converge.  The  old, 
familiar,  and  perhaps  inefficient  interval -halving  method  was  used  for 
this  procedure.  Further  discussion  concerning  this  procedure  is  given 
in  Chapter  VI. 
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V  Results  —  Flow  Over  a  Circular  Cylinder 

Finite  element  solutions  are  presented  for  flow  over  a  unit 
circular  cylinder  for  Mach  numbers  within  the  range  of  Oi  -  Mc t 

These  solutions  are  compared  with  exact  solutions  for  incompressible 
flow  and  with  solutions  obtained  from  other  approximate  methods  for 
compressible  flow.  The  problem  of  incompressible  flow  without  circu¬ 
lation  was  solved  using  three  different  elements  representing  one 
nonconforming  and  two  conforming  approximations.  The  two  conforming 
elements  were  used  to  solve  the  incompressible  problem  with  circula¬ 
tion.  The  compressible  problem  was  solved  using  the  new  conforming 
element. 


Incompressible  Flow  Without  Circulation 

Geometry.  The  entire  flowfield  was  discretized  as  shown  in 
Fig  5  for  a  typical  set  of  discretization  parameters.  Only  a  quarter 
of  the  field  is  necessary  because  of  syometry,  but  since  the  addition 
of  circulation  requires  the  entire  field,  then  it  was  discretized  for 
all  of  the  incompressible  problems.  A  computer  subroutine  was  written 
to  automatically  discretize  the  field  from  the  parameters  (  *FF,  Nr. 
Ns).  Rpp  is  the  radius  of  the  farfield  boundary,  the  number  of 

element  rings,  and  the  number  of  angular  sectors  in  each  ring 

(see  Fig  5) . 

The  location  of  the  outer  boundary  depends  upon  the  boundary 
condition  imposed  there.  When  the  farfield  expression  given  by  eq  40 
is  used,  can  be  as  close  as  3  radii  from  the  cylinder.  If  other 

boundary  conditions  are  used,  for  example  O  *  then  ftpp  has 
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to  be  extended  beyond  10  radii  before  the  velocity  profile  on  the 
cylinder  surface  is  unaffected.  The  boundary  condition  *  O  does 

not  produce  adequate  velocity  distributions  when  far field  boundaries 
are  located  within  10  radii  of  the  surface.  This  behavior  was  expected 
and  affects  the  velocity  distribution  significantly  in  the  region  of 
maximum  velocity.  For  this  reason,  and  to  preclude  the  necessity  of 
considering  large  flow  domains ,  the  farf ield  boundary  condition  given 
by  eq  40  is  used  exclusively. 

Comparison  of  Three  Sector  Elements.  Before  examining  the  effect 
of  element  refinement  on  convergence,  a  comparison  among  the  three 
finite  element  approximations  and  the  exact  solution  is  given.  The 
solutions  for  a  coarse  flowfield  discretization  of  only  12  elements 
(3  rings  and  4  sectors)  a^e  shown  in  Fig  7  for  the  three  trial  functions 
described  in  Chapter  III.  The  nonconforming  approximation  appears  to 
be  better  than  either  of  the  other  two  conforming  ones  for  the  dis¬ 
cretization  used.  The  noncomforming  approximation  obtained  from  this 
coarse  discretization  actually  compares  well  with  solutions  obtained 
from  the  other  two  conforming  elements  when  a  more  refined  discreti¬ 
zation  is  used.  One  might  conclude,  on  the  basis  of  these  results 
and  observations,  that  sector  element  (1)  is  the  most  desirable 
element;  however,  this  choice  is  a  bad  one.  For  the  two  conforming 
elements,  the  approximations  improve  when  the  discretization  is  re¬ 
fined,  but  they  do  not  improve  for  the  nonconforming  element.  In 
fact,  the  error  in  the  approximations  becomes  greater  when  the  dis¬ 
cretization  is  refined. 

One  reason  for  the  decline  in  accuracy  for  the  nonconforming 
element  is  due  to  the  error  caused  by  not  satisfying  the  required 
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continuity  of  potential  across  inter-element  boundaries.  For  coarse 
discretizations,  which  mean  relatively  few  such  boundary  lines,  the 
error  is  comparatively  small.  However,  as  the  discretization  is 
refined  additional  error  occurs  and  the  solution  degenerates.  Further 
refinement  compounds  the  problem  which  leads  to  the  rejection  of  sector 
element  (1) ,  except  for  the  most  coarse  discretizations .  There  are 
formulation  adjustments  which  can  be  made  to  improve  the  nonconforming 
solutions  significantly  for  finer  discretizations.  However ,  the 
resulting  solutions  are  not  as  good  as  those  obtained  from  the  con¬ 
forming  elements  and  these  adjustments  will  not  be  discussed. 

Discretization  Effects.  Tangential,  velocity  distributions 
v6(v=i,e)  obtained  from  each  conforming  element  are  compared  with 
the  exact  solution  as  shown  in  Figs  8  and  9 .  The  four  finite  element 
solutions  in  each  figure  are  for  four  different  discretizations  of 
p  .  Additionally,  these  discretizations  represent  refinement  of 
the  element  angular  size  parameter  ^  for  a  fixed  value  of  the  element 
radial  size  parameter  .  These  solutions  are  presented  to  show  the 
difference  in  convergence  properties  between  the  two  conforming 
elements  as  a  function  of  discretization  refinement.  By  comparing  the 
solutions  in  Figs  8  and  9  with  the  exact  solution  it  appears  that 
refinement  of  alone  is  sufficient  to  achieve  convergence  with  element 
(3)  ,  the  new  element.  It  also  appears  that  refinement  of  (2  alone  is 
not  sufficient  to  achieve  convergence  with  element  (2),  the  bilinear 
element.  For  each  discretization  the  velocity  distribution  obtained 
with  element  (3)  is  closer  to  the  exact  distribution  than  the  distri¬ 
bution  obtained  from  element  (2).  These  differences  are  due  to  the 
diverse  nature  of  the  trial  function  for  each  element. 
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Figure  9  -  Convergence  of  Tangential  Velocity  for  Sector 
Element  (3)  for  Refinement  of  Q  Alone 
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To  further  illustrate  the  difference  in  convergence  properties 
between  the  two  conforming  elements,  consider  two  variations  in  the 
discretizations  of  siF.  First,  the  effect  of  refining  only  the 
radial  parameter  o(  is  examined.  The  number  of  angular  sectors  and 
location  of  the  farfield  boundary  remain  fixed.  Table  I  shows  the 
effect  on  the  solution  of  the  tangential  velocity  at  0  :  t  TT /  2. 
for  refinement  of  o(  .  The  points  Q  z  tlT/a.  are  selected  because 
the  point -wise  error  there  is  a  maximum  as  shown  in  Figs  8  and  9. 

Table  I  -  Point -wise  Error  for  Conforming 

Sector  Elements  with  Radial  Refinement* 


Number 

of 

Rings 

1 

Outer  Radii  of  Ring 

Ve(l,TT/2.) 

Per  Cent  Error 

1  2  3  4  5  6 

Element  (3) 

Element  (2) 

3 

2.0  4.0  7.8 

0.95 

4.46 

4 

1.4  2.0  4.0  7.8 

0.95 

2.49 

5 

1.4  2.0  2.8  4.0  7.8 

0.95 

1.99 

6 

1.4  2.0  2.8  4.0  5.6  7.8 

0.95 

1.87 

_  .  i 

4 

2.0  2.8  4.0  7.8 

0.95 

-1 

3.98 

4 

2.0  4.0  5.6  7.8 

0.95 

4.35 

1 

*  Number  of  Angular  sectors  equals  16. 

The  percent  error  is  based  on  the  exact  solution  and  is  presented  for 
both  conforming  elements.  First,  note  that  the  error  in  velocity 
V9(i,tr/i)  for  sector  element  (3)  is  constant  at  0. 95%  regardless  of 
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the  number  of  rings  or  the  relative  size  arrangement  of  those  rings. 
Thus,  radial  refinement  is  not  required,  and  does  not  improve  the 
results  obtained  with  this  element.  The  reason  for  this  behavior  is 
due  to  the  functional  dependence  of  f  in  the  trial  function  Since 
the  exact  solution  behaves  as  CoZd/v*  *  then  r^e  trial  function 
contains  the  correct  form  for  the  radial  variable.  This  behavior 
suggests  interesting  possibilities  for  other  types  of  problems  which 
will  be  discussed  in  more  detail  in  a  later  section  on  recormendation 
for  further  work.  Secondly,  note  that  refinement  of  parameter  oi 
for  sector  element  (2)  significantly  affects  the  solution  For  this 
element  the  radial  dependence  of  y*  in  the  trial  function  is  different 
from  the  exact  form;  therefore,  to  achieve  convergence,  refinement  of 
is  necessary. 

Three  other  interesting  observations  can  be  made  about  element  (2) 
from  the  data  in  Table  I.  First  of  all,  radial  refinement  alone  will 
not  drive  the  error  to  zero.  For  a  given  number  of  angular  divas ions 
rhere  is  a  certain  number  of  rings  beyond  which  the  solution  does  not 
get  closer  to  the  exact  one  This  behavior  is  detected  by  observing 
the  rate  at  which  the  error  decreases  as  more  rings  are  added. 

Secondly,  if  only  one  additional  ring  is  to  be  added  to  some  base-line 
discretization  (represented  by  3  rings  in  Table  I)  then  the  location 
of  that  ring  is  crucial.  For  example,  approximately  447  reduction  in 
error  is  achieved  bv  dividing  the  first  ring  into  two  rings  of  pro¬ 
portional  thickness  (as  illustrated  by  line  two  in  the  table) .  Similar 
divisions  of  the  second  and  third  rings  result  in  approximately  11% 
and  2%  reduction,  respectively  (as  illustrated  by  the  last  two  lines 
in  the  table).  Thirdly,  for  all  of  the  discretizations  shown,  the 
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error  obtained  from  the  new  element,  element  (3),  is  significantly 
less  than  from  the  bilinear  element .  element  (2) . 

Next,  the  effect  of  refining  only  the  angular  size  parameter 
is  examined.  Refinement  of  this  parameter  is  required  of  both  con¬ 
forming  elements  to  achieve  convergence.  The  number  of  rings  and 
location  of  the  far field  boundary  remain  fixed.  The  tangential  velo¬ 
city  at  ±TT/z  is  presented  in  Table  II  for  different  angular 
divisions . 

Table  II  -  Point -wise  Error  for  Conformal 

Sector  Elements  with  Angular  Refinement** 


Number 

Ue<7,  tir/z) 

Percent  Error 

of  Sectors 

Element  (3) 

Element  (2) 

Element  (3) 

Element  (2) 

4 

1  7038 

1.6484 

14.81 

17.58 

8 

1 . 9244 

1.8569 

3.78 

7.16 

16 

1.9810 

1.9108 

.95 

4.46 

32 

1.9952 

1.9244 

.24 

3.78 

48 

- 

1.9979 

1.9269 

.11 

3.66 

**Nimber  of  rings  equals  3. 


The  error  in  velocity  at  ©  =  -^2  for  the  new  element  can  be  shown  to 
decrease  by  0(t)  The  bilinear  element  converges  rapidly  for  a 
couple  of  refinements  in  ^  ,  but  then  a  point  is  reached  where  further 
refinement  produces  very  little  change  of  results.  This  same 
behavior  was  observed  in  Table  I  for  radial  refinement.  Thus,  for 
element  (2)  both  the  radial  and  angular  size  parameters  must  be 


refined  simultaneously  with  some  other  parameter  such  as  element 
aspect  ratio  held  constant. 

Incompressible  Flow  With  Circulation 

Most  of  the  discretization  questions  discussed  for  flow  without 
circulation  also  apply  when  circulation  is  present.  The  required  far- 
field  boundary  location,  the  effect  of  element  refinement,  and  the 
characteristic  differences  between  the  two  conforming  elements  are 
unchanged  by  the  addition  of  circulation.  The  only  remaining  dis¬ 
cretization  question  that  needs  to  be  answered  is  whether  the  relative 
arrangement  of  elements  about  the  stagnation  point  has  an  impact  on 
the  value  of  circulation  predicted. 

Circulation.  In  general,  the  value  of  circulation  obtained  from 
the  solution  depends  upon  the  relative  location  of  the  stagnation  point 
within  the  element  containing  it.  To  quantify  this  statement  Table  III 
presents  the  error  in  predicted  circulation  for  a  number  of  discreti¬ 
zations.  The  data  in  this  table  was  obtained  with  the  use  of  element 
(3)  with  the  farfield  boundary  located  at  ftfF-S  and  stagnation  point 
at  ©=-TT/6  ■ 
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Table  III  -  Error  in  Circulation  for  Sector  Element  (3) 
for  Stagnation  Point  Q  r  -"TT/6 


Angular 

Sectors 

Number 

Element 

Number 

Nodes 

"  ’ 

%  Error  in 
Circulation 

8 

16 

27 

11.6 

12 

24 

30 

24.0 

16 

32 

51 

5.5 

20 

40 

63 

4.6  i 

24 

48 

75 

11.6 

32 

6^ 

99 

2.9 

40 

80 

123 

2.2 

The  error  in  circulation  appears  to  behave  irratically  as  the  number  of 
angular  sectors  is  increased.  One  would  expect  the  error  to  decrease 
as  the  parameter  ^  is  refined.  The  reason  for  the  erratic  behavior 
is  due  to  the  relative  difference  in  locations  of  the  stagnation  point 
within  the  appropriate  stagnation  element  as  ^  is  refined.  When 
location  differences  are  accounted  for,  then  the  error  in  circulation 
does  decrease  as  ^  is  refined.  For  example,  for  angular  divisions  of 
12  and  24  the  stagnation  point  lies  on  a  node.  Doubling  the  number  of 
angular  divisions  from  12  to  24  reduces  the  error  in  circulation  by 
approximately  50%;  although,  the  magnitude  of  error  is  greater  than  for 
solutions  with  fewer  angular  divisions  (i.e.  larger  elements).  Also, 
doubling  the  number  of  angular  divisions  from  8  to  16  and  again  to  32 
shows  an  error  reduction  of  about  50%  for  each  doubling.  The  same 
observation  holds  true  for  the  division  from  20  to  40  sectors.  So, 
one  may  wonder,  what  causes  the  erratic  behavior? 
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A  search  was  made  to  determine  if  an  optimum  point  exists  within 
the  stagnation  element  such  that  the  error  in  circulation  is  minimal. 


Such  a  point  was  found  to  be  the  angular  centroid  of  the  element..  Once 
the  stagnation  point  is  established,  then  the  flowfield  should  be  dis¬ 
cretized  to  permit  the  stagnation  point  to  be  located  at  the  angular 
centroid  of  the  stagnation  element.  For  instance,  in  order  for  a 
stagnation  point  at  £c~TT/£  to  occur  at  the  proper  location  in  the 
element,  either  6,  18,  or  30  angular  divisions  should  be  used.  Table 
IV  gives  the  results  from  both  conforming  elements  using  these  angular 
divisions  and  4  rings  of  elements. 


Table  IV  -  Error  in  Circulation  for  Conforming 
Sector  Elements  with  Stagnation 
Point  Tr/to  Located  at  the  Optimum 
Point. 


Angular 

Sectors 

Number 

Nodes 

70  Error  in  Circulation 

Element  (3) 

Element  (2) 

6 

28 

0.0144 

3.797 

18 

95 

0.0080 

0.597 

30 

155 

0.0032 

0.597 

It  is  readily  apparent  by  comparing  the  results  presented  in  Tables 
III  and  IV,  that  the  error  obtained  for  circulation  strongly  depends 
on  where  the  stagnation  point  is  located  in  the  stagnation  element.  The 
error  is  maximum  when  the  stagnation  point  coincides  with  a  node,  and 
it  is  mimimum  when  the  point  lies  half  way  between  two  nodes.  It  is 
also  apparent  that  for  a  given  discretization  the  new  element  provides 
by  far  the  best  approximation,  although  acceptable  results  are  obtained 


from  the  bilinear  element  for  the  coarse  discretizations  presented. 

Velocity  Distribution.  Figures  10  and  11  give  velocity  distri¬ 
butions  obtained  fran  each  of  the  conforming  elements  for  the  discreti¬ 
zations  presented  in  Table  IV.  The  exact  solutions  are  plotted  for 
comparison  purposes.  Accurate  estimates  of  the  velocity  profile  are 
achieved  whenever  the  circulation  can  be  predicted  correctly. 

Compressible  Flow 

Iterative  Algorithm.  The  iterative  solution  algorithm  for  con¬ 
stant  coefficients  defined  by  eq  44  converged  for  subsonic  flow  and 
even  for  flows  where  the  local  Mach  number  exceeded  one.  It  is  believed 
that  a  minor  disagreement  between  the  onset  of  transonic  flow  and  the 
simultaneous  occurrence  of  a  change  in  the  differential  equation  from 
elliptic  to  hyperbolic,  accounts  for  the  convergence  for  slightly 
supercritical  flow.  Since  the  small-disturbance  equation  is  not  totally 
applicable  for  compressible  flow  over  a  cylinder,  and  since  it  is 
locally  linearized  by  constant  local  coefficients,  then  it  does  not 
change  types  at  exactly  the  same  for  which  the  local  Mach  number 

exceeds  unity.  The  local  Mach  number  must  exceed  unity  to  the  extent 
that  an  element  becomes  supersonic  (engulfed  within  the  supersonic 
bubble)  before  the  equation  changes  types.  The  occurrence  was  created 
by  formulation  assumptions  and  would  not  occur  if  the  nonlinearity  were 
treated  in  a  more  exact  manner.  When  the  differential  equation  for  one 
or  more  elements  becomes  hyperbolic  the  iterative  solution  diverges . 

A  few  reported  "upwinding"  algorithms  were  attempted,  but  none  proved 
successful  and  the  attempt  was  abandoned  for  the  cvlinder.  The  new 
"upwinding"  technique  described  for  the  airfoil  was  not  tried  and  mav 
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Figure  10!)  -  Perturbation  Velocity  Distribution  from  Sector 
Element  (3)  for  Discretization  (2) ,  Table  IV 
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prove  useful  for  transonic  flows. 

Critical  Mach  Number.  The  critical  Mach  number  =  0.403  was 
predicted  for  a  discretization  of  64  elements  with  80  nodes.  This 
value  compares  well  with  that  obtained  by  Imai  (Ref  30) ,  who  used  a 
Janzen-Rayleig.h  method  of  third  order.  Table  V  shows  a  comparison  of 
critical  Mach  numbers  reported  by  several  investigators.  Greenspan 
(Ref  32)  used  finite  difference  methods  in  conjunction  with  a  varia¬ 
tional  principle.  Habashi  (Ref  31)  used  linear  triangular  elements  as 
discussed  in  Chapter  I. 

Table  V  -  Comparison  Critical  Mach  Numbers 


Source  of  Results 

%  Difference* 

This  Report 

0.403 

-0.32 

Imai  (Ref  30) 

— 

3rd-0rder 

0.4043 

0.00 

2nd-Order 

0.4090 

1.16 

Ist-Order 

0.4206 

4.03 

Habashi  (Ref  31)  ^ 

0.40-0.42 

-1.06  to  +3.88 

Greenspan  (Ref  32)  j 

0.404 

-0.07 

*  %  Difference  based  upon  Imai's  3rd-order  approximation. 


Velocity  Distributions .  The  velocity  distributions  for  M#  =0.3 
and  0.4  are  presented  in  Fig  12.  Imai’s  third-order  acciirate  results 
are  presented  for  comparison  purposes.  Agreement  is  good  over 
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approximately  60  degrees  of  a  quadrant.  Differences  are  more  pro¬ 
nounced  in  the  neighborhood  ofTT/2.  It  is  in  this  region  that  finite 
element  results  for  smaller  elements  are  not  as  good  as  those  reported. 
It  is  believed  that  the  error  resulting  from  using  the  small-disturbance 
equation  is  responsible  for  the  deterioration  when  smaller  elements  are 
used.  Velocity  perturbations  in  this  region  approach  or  exceed  . 


j 


©  (radians) 

Tangential  Velocity  Distribution  on  the  Cylinder  Surface  Without  Circulation 
(Compressible  Flow) 
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VI  Results  of  Flow  Over  a  Nonlifting  Airfoil 

Finite  element  results  are  presented  for  a  symmetric,  6%- thick, 
parabolic-arc  dirfoil.  Pressure  coefficient  distributions  on  the 
upper  surface  of  the  airfoil  are  computed  from  the  finite  element 
solutions  of  the  potential  function  in  domain  Slf.  Appendix  D  con¬ 
tains  the  computational  details  of  this  process.  Distributions  of 
the  pressure  coefficient  are  presented  for  a  range  of  Mach  numbers 
from  zero  up  through  the  transonic,  mixed- flow  regime.  Where  possible, 
these  distributions  are  compared  with  experimental  data  and  distribu¬ 
tions  computed  from  classical  thin-airfoil  theory,  finite  difference 
methods,  and  other  finite  element  methods.  Convergence  properties  of 
the  iterative  solution  algorithm  are  investigated.  In  addition,  con¬ 
vergence  properties  of  the  pressure  distribution  as  a  function  of 
discretization  refinement  and  location  of  the  farfield  boundary  are 
established. 

Discretization  Effects 

When  the  flow  domain  is  discretized  several  decisions  must 
be  made  concerning  the  flowfield  size  and  arrangement  of  elements.  For 
instance,  one  must  decide  where  to  place  the  farfield  boundary.  Should 
it  be  relatively  close  to  the  airfoil  or  far  from  it?  How  many 
elements  should  be  used?  If  smaller  elements  are  needed  to  improve 
the  accuracy  of  desired  results,  then  how  should  the  discretization 
be  refined?  Should  the  discretization  of  the  entire  flow  domain  be 
refined?  Will  refinement  in  a  subregion  of  £Lp  improve  desired 
results?  If  so,  then  which  subregion?  Answers  to  these  and  related 
questions  are  given 
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Far field  Boundary  Location.  The  far field  boundary  location  of 
flow  domain  £lp  can  be  characterized  by  two  parameters  X^x  and 
Ymax  as  shown  in  Fig  13  The  solution  for  the  airfoil  problem  as 
formulated  could  be  obtained  for  any  combination  of  X^*  and  Ymx. 

One  could  select  relatively  large  values  for  these  parameters,  but 
from  a  computational  viewpoint  this  approach  may  be  prohibitive.  One 
would  like  to  select  the  smallest  possible  region  to  reduce  the  cost 
of  computing  a  solution  without  sacrificing  desired  accuracy.  Thus, 
lower  permissible  bounds  of  Xma*  and  Y^*  were  estimated  from  a  series 
of  solutions  obtained  from  elements  of  fixed  size.  First,  a  baseline 
solution  was  computed  for  a  given  value  of  Xmax  and  Y^ax  By  adding 
additional  elements  to  the  outer  perimeter  of  the  flow  domain  a 
succession  of  solutions  was  obtained  for  increasing  values  of  X^ax 
and  Y roQx.  By  holding  the  element  size  constant  and  varying  only  the 
number  of  elements  used,  the  effect  on  the  pressure  distribution  as 
a  function  of  only  the  location  of  the  farfield  boundary  could  be  as¬ 
sessed.  This  process  was  done  twice,  once  to  get  a  bound  on  X^ax  for 
Ymax  held  constant  and  vice  versa.  Figures  14  and  15  show  the  effect 
on  the  pressure  distributions  obtained  by  this  procedure  for  Mco  ~  0* 
There  is  nothing  magic  about  selecting  the  number  3/2  for  the  parameter 
held  fixed.  Any  greater  value  could  be  chosen  to  demonstrate  trends. 
The  table  in  Figs  14  and  15  gives  the  numerical  values  of  the  pressure 
coefficient  at  the  midchord  of  the  airfoil.  From  the  trend  of  this 
data  the  pressure  coefficient  converges  from  below  as  X^ax  and  Ymax 
are  increased.  This  means  that  the  perturbations  velocities  ^con¬ 
verge  from  above  as  the  farfield  boundary  is  extended.  From  the  data 
presented  in  these  two  figures  there  appears  to  be  a  farfield  boundary 
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location  beyond  which  pressure  distributions  are  unaffected  as  the 
boundaries  are  extended. 

Other  descriptions  of  convergence  are  provided  in  Figs  16  to  19 
which  include  the  effect  of  Mach  number.  The  pressure  coefficient 
in  these  figures  is  normalized  with  respect  to  its  converged  value. 
Figures  16  and  17  are  for  the  1/4-chord  point  (also  3/4-chord)  while 
Figs  18  and  19  are  for  the  midchord  point.  Convergence  trends  are 
the  same  at  these  three  points  and  as  demonstrated  in  Figs  14  and  15 
for  incompressible  flow,  the  trend  is  similar  for  the  entire  distribu¬ 
tion. 

What  conclusions  can  be  made  concerning  the  farfield  boundary 
location  from  these  figures?  First,  the  curves  show  that  farfield 
boundaries  more  than  3  chord  lengths  from  the  origin  will  not  appreci¬ 
ably  improve  results.  In  fact,  farfield  boundaries  could  be  chosen 
closer  with  little  effect.  For  example,  less  than  4%  error  occurs 
when  farfield  boundaries  are  located  as  close  as  3/2  chord  lengths. 
Secondly,  the  effect  of  Mach  number  on  the  lower  bound  for  Xmax  is 
opposite  to  that  for  Yn]ax.  For  a  desired  degree  of  accuracy  Xmax 
can  be  taken  smaller  while  Yrriax  should  be  larger  as  Mach  number  in¬ 
creases.  This  behavior  is  not  surprising  in  light  of  the  Prandtl- 
Glauert  transformation  which  stretches  the  y  -variable  as  Mach  number 
increases.  Thus,  perturbations  need  more  distance  in  the  ^  -direction 
to  dissipate. 

Farfield  Boundary  Conditions .  Associated  with  the  farfield 
boundary  location  is  the  condition  imposed  there.  For  nonlifting  flow 
the  farfield  boundary  condition  for  the  potential  is  given  by  eq  50 
with  r  r  O  .  The  first  and  last  terms  of  that  equation  can  be  expressed 
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as 


4>Ff  =  4>pf  +•  4>ff  (( 

inl  ju 

where  (pff  is  of  higher  order  than  (p^p  ,  as  previously  discussed.  Now? 
suppose  that  the  above  boundary  condition  is  modified  by  a  constant 


and  written  as 


<!fF  »  K  (  4>pf  +  4>FF  ) 


When  (<  -  1  ,  the  condition  is  identical  to  that  of  Klunker,  eq  50 

For  any  other  value  of  the  doublet  strength  expressed  in  eq  50  is 

j  NL 

proportionally  modified.  Recall  that  the  term  i comes  from  an 
integration  of  over  the  entire  domain  which  means  that  it  depends 

upon  the  solution.  When  this  term  is  kept  in  the  far field  boundarv 
condition  to  compute  the  solution  for  the  n-th  iteration  of  (£)  ,  then 
velocity  from  the  previous  iteration  is  used . 

Now,  consider  the  effect  on  the  pressure  coefficient  as  K  is 
varied  from  0  to  1,  as  shown  in  Fig  20  for  incompressible  flow’  K  ~  o 
is  a  bounding  condition  for  the  distribution  of  pressure  coefficient 
Cp  Increasing  K  has  the  effect  of  translating  the  CLp  distribu¬ 
tion  uniformly  toward  more  negative  values.  Since  convergence  of 
Cp  is  from  below,  as  will  no  shown  in  the  next  section,  then  in¬ 
creasing  K  can  produce  very  accurate  results  for  relatively  coarse 
discretizations  (i.e.  largo  elements).  In  fact,  for  any  real  is?  i.e 
discretization  that  might  be  used  a  value  of  K  can  he  found  that  will 
provide  C_p  approximations  that  accurately  match  known  or  experi¬ 
mental  results.  All  of  the  results  presented  in  this  chapter  are  for 
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For  compressible  flow  the  effect  of  keeping  the  nonlinear  far- 
field  boundary  term  is  not  significant  (Ref  13) .  Since  this  term  is 
of  higher  order  (Ref  18)  its  contribution  is  negligible  particularly 
for  subsonic  flow  and  even  for  Mach  numbers  which  produce  transonic 
flow.  For  example,  the  average  change  in  the  solution  for  the  pressure 
coefficient  with  the  nonlinear  boundary  term  included  is  only:  0 . 097, 
for  0.5,  0.327-  for  Meo=  0.7,  and  0.85%  for  0.8.  Due  to 

the  computational  nature  of  (ppf  and  its  minimal  effect  on  the  solution, 
it  may  be  neglected  in  favor  of  reducing  computer  time. 

No  other  type  of  farfield  boundary  conditions  were  studied;  how¬ 
ever,  several  other  possibilities  could  be  tried  and  will  undoubtedly 
be  the  basis  for  further  study. 

Element  Refinement.  In  order  to  study  the  point -wise  convergence 
of  the  pressure  coefficient  as  element  size  decreases,  all  other  dis¬ 
cretization  and  solution  parameters  were  held  fixed  while  element  size 
was  varied.  Figure  13  shows  a  typical  discretization  of  the  flow- 
field  governed  by  parameters  NDX,  NDY,  and  NDXA  which  are  defined  in 
the  figure.  The  farfield  boundary  was  located  at  Xmax  =  Ymax  =1.50.  • 
Element  refinement  was  done  4  ways: 

(1)  All  elements  in  the  flowfield  were  uniformly  decreased 
in  size. 

(2)  Only  elements  above  the  airfoil  were  decreased  in  width 
only  while  those  for  I  %\*oS  were  held  fixed  (see  Fig  13) . 

(3)  Elements  for  |X|»*.r  were  refined  in  width  only  while 
all  other  discretization  parameters  were  held  fixed 

(i.e.  changing  AX  for  /*/  ~>0S  while  for 


83 


and  for  y  >  C>  are  fixed ;  see  Fig  13)  . 

(4)  All  elements  were  refined  in  height  (i.e.  giving  more 
layers  of  elements)  while  all  other  discretization 
parameters  were  held  fixed. 

Figure  21  shows  the  effect  of  uniform  element  refinement  (i.e.  proce¬ 
dure  1  above)  on  the  solution  for  the  pressure  coefficient  for  tAeo'O. 
The  exact  curve  computed  from  thin-airfoil  theory  is  also  shown  for 
comparison  purposes. 

IVro  observations  can  be  made  from  Fig  21.  First,  the  step  function 
for  each  of  the  approximations  intersects  the  exact  curve  at  some  point 
in  the  interval  of  the  step  size.  Thus,  even  the  approximation  from 
the  most  coarse  discretization  is  not  unrealistic.  Secondly,  the 
point-wise  error  over  most  of  the  chord  decreases  with  a  corresponding 
decrease  in  element  size.  In  an  average  sense,  the  error  in  pressure 
decreases  with  element  refinement.  Further,  the  point-wise  convergence 
at  the  midchord  appears  to  be  very  rapid  compared  to  that  at  the 
leading  and  trailing  edges.  The  description  of  pressure  coefficient 
at  the  leading  and  trailing  edges  improves  as  the  elements  become 
smaller,  but  improvement  is  slower  than  for  other  points  on  the  air¬ 
foil.  Special  finite  element  treatment  may  be  needed  for  regions  near 
the  singular  points  to  improve  the  convergence  there.  The  results  in 
Fig  21  show  that  the  solution  trend  is  correct,  in  that,  it  approaches 
the  exact  solution  as  element  size  approaches  zero. 

Refinement  of  elements  only  in  the  region  over  the  airfoil 
profile  (i.e.  procedure  2)  has  a  slightly  different  effect  on  the 
pressure  distribution  than  does  uniform  refinement  of  all  elements 
in  the  flowfield.  The  coefficients  of  pressure  for  refinement  of 
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disc?  atization  parameter  NDXA  alone  are  presented  in  Fig  22.  Note 
that  the  actual  shape  of  the  pressure  distribution  is  more  well  de¬ 
fined  with  smaller  elements  due  to  the  step-function  nature  of 
pressure,  but  the  pressure  at  the  mid-point  of  each  step  does  not 
show  substantial  deviation  from  the  exact  solution  with  refinement. 

The  greatest  change  with  refinement  of  NEKA  occurs  near  the  leading 
and  trailing  edges  where  better  definition  of  the  pressure  coefficient 
is  obtained  as  the  elements  become  smaller.  The  change  near  the  peak 
is  almost  insignificant  by  comparison. 

The  difference  in  the  behavior  of  pressure  for  the  two  refinement 
procedures  described  thus  far  suggests  a  discretization  technique. 
First ,  discretize  -ft.  p  with  relatively  large  elements  everywhere  and 
refine  the  element  size  until  there  is  little  change  in  the  pressure 
coefficient  at  the  midchord.  Then,  for  better  definition  of  the  dis¬ 
tribution  refine  only  the  NEKA  parameter.  Perhaps  an  improvement  on 
this  technique  would  be  the  use  of  variable  size  elements  with  the 
smaller  elements  located  near  the  leading  and  trailing  edges. 

Refining  elements  beyond  the  airf oil  leading  and  trailing  edges 
(i.e.  procedure  3,  refining  parameter  NDX  only)  has  little  effect  on 
the  pressure  distribution  except  in  the  two  elements  containing  the 
leading  and  trailing  edges  respectively;  refer  to  Fig  23.  The  varia¬ 
tion  of  pressure  that  occurs  in  these  two  elements  does  not  appear  to 
be  directly  related  to  the  change  made  in  the  refined  elements.  The 
fact  that  dip  changes  at  all  is  due  to  the  decrease  of  size  in 
elements  adjacent  to  the  leading  and  trailing  edge  elements.  It 
clearly  points  out  that  "for  points  farthest  removed  from  the  elements 
refined,  the  smaller  the  change  the  refinement  makes."  This  "truism" 
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is  detected  in  Fig  23  by  noting  the  decreasing  rate  of  change  in 
pressure  as  the  tnidchord  is  approached  from  either  direction. 

Refinement  of  NDY  with  all  other  discretization  parameters  held 
fixed  (i.e.  procedure  4)  significantly  affects  the  pressure  distribu¬ 
tion  as  shown  in  Fig  24.  Changes  occur  rapidly  near  the  peak,  but 
are  slower  near  the  leading  and  trailing  edges.  The  width  of  the 
elements  as  well  as  the  height  needed  to  be  refined  near  these  points 
to  improve  the  resulting  approximation. 

Incompressible  Results 

Figure  25  shows  a  comparison  between  the  finite  element  solution 
for  the  pressure  coefficient  and  the  exact  curve  computed  from  thin 
airfoil  theory.  The  pressure  coefficient  is  a  step  function,  but  when 
the  midpoint  value  in  each  step  interval  is  plotted  the  comparison 
with  the  exact  solution  shows  good  agreement.  The  finite  element 
solution  was  obtained  vising  discretization  parameters.-  Xmax  =  1.5, 
Ymax  =  2.0,  NDX  =  8,  NDXA  =  18,  NDY  =  8. 

Another  interesting  comparison  is  shown  in  Fig  26  where  the 
finite  element  solution  is  compared  with  a  finite  difference  solution 
obtained  from  Ref  (40) .  The  finite  difference  solution  was  obtained 
using  a  constant  grid  step  size  AX  =  A'd  =  0.125  everywhere.  The 
finite  element  results  were  obtained  using  the  same  step  size  over  the 
airfoil  surface  (i.e.  AX^  =  0.125),  but  the  other  element  dimensions 
were  twice  as  large  as  that  used  in  the  finite  difference  calculations 
(i.e.  AX  =  0.25,  A'j)  =  0.25).  The  degrees -of- freedom  for  the  two 
solutions  were  119  nodes  for  finite  element  and  350  grid  points  for 
finite  difference.  The  accuracy  achieved  by  the  finite  element  method 
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is  clearly  better  than  for  finite  difference. 

Compressible  Subsonic  Results 

The  flow  is  subsonic  at  all  points  in  the  flow  domain  Si-  p 
as  long  as  the  coefficient  of  $xX  •*-n  the  governing  differential 
equation  remains  positive.  When  the  free-stream  Mach  number  is  in¬ 
creased  slightly  above  a  critical  value,  then  the  coefficient  of 
becomes  negative  in  a  small  subregion  of  £Lf  and  transonic  flow  begins. 
Before  discussing  results  for  compressible  subsonic  flow,  the  conver¬ 
gence  behavior  of  the  iterative  solution  algorithm  for  Mach  numbers 
through  the  subsonic  regime  up  to  the  onset  of  transonic  flow  is  con¬ 
sidered. 

Convergence  of  Iterative  Scheme.  Recall  that  the  approximate 
solution  of  the  governing  differential  equation  was  written  in 
iterative  form  as  expressed  by  eq  58.  Also  recall  that  convergence  of 

the  iterative  solution  scheme  was  governed  by  the  criteria  specified 

.  -V 

in  eq  45.  For  £  r  0,5  XlO  the  number  of  iterations  required  to 
achieve  convergence  of  eq  58  is  presented  in  Figs  27  and  28  as  a 
function  of  Mach  number.  For  the  data  in  these  figures,  solutions 
were  started  by  setting  all  nodal  values  of  potential  equal  to  zero. 
For  low  subsonic  flow  (p.Of:  Moo  -  0,5])  the  number  of  iterations  re¬ 
quired  for  convergence  of  the  iterative  scheme  is  five  or  less. 

Fewer  iterations  were  required  when  solutions  were  started  from  better 
initial  guesses.  For  example,  fewer  iterations  than  given  in  Fig  27 
were  required  to  obtain  a  solution  for  the  iterative 

scheme  was  started  from  the  solution  for  Moo-  0.  4  •  As  the  Mach 
number  was  increased  to  the  slope  of  the  convergence  curve 

increased,  but  convergence  still  occurred  in  seven  or  less  iterations. 
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At  approximately  M®  =  0.85  a  sharp  increase  in  the  number  of 
iterations  occurred  as  indicated  by  the  vertical  segment  of  the  curve 
in  Fig  27.  At  this  Mach  number  a  small  supersonic  bubble  has  formed 
in  the  flow  above  the  midchord  of  the  airfoil  and  the  equation  has 
locally  changed  type.  For  0<Q5  the  iterative  solution  scheme 

no  longer  converged  in  the  sense  of  eq  45  for  small  values  of  £  . 

What  happened  to  the  solution  of  the  potential  function  for  each 
iteration  is  illustrated  in  Fig  28.  In  this  figure  the  percent  change 


in  the  nodal  value  of  the  potential  function  at  the  midchord  was  com¬ 
puted  frcm  and  plotted  as  a  function  of  the 

iteration  number  for  Mach  numbers  frcm  0.50  to  0.86.  The  solution  was 


initiated  by  setting  the  potential  function  equal  to  zero  at  all  nodes. 
As  observed  in  Fig  27,  for  low  subsonic  flow  the  potential  converges 
rapidly,  but  for  the  potential  function  does  not  converge 

at  all.  Initially  the  solution  gives  the  appearance  of  converging  as 
it  does  for  lower  Mach  numbers ,  but  then  an  iteration  is  reached 
(iteration  6)  where  the  apparent  convergence  trend  begins  to  reverse 
itself.  This  behavior  occurs  because  the  governing  equation  has  locally 
changed  type  from  elliptic  to  hyperbolic  in  a  small  region  above  the 
michord.  Further  iterations  produced  what  appears  to  be  a  diverging 
solution  scheme.  In  actuality  the  solution  does  not  diverge,  but 
cycles  back  through  the  "bucket"  shape  shown  in  Fig  28.  The  behavior 
for  larger  Mach  numbers  was  similar  to  that  shown  for  0.85-  In 

general,  the  "bucket"  shape  moves  to  the  left  and  upward  indicating 
that  fewer  iterations  are  required  before  the  solution  develops  to 
the  stage  where  transonic  flow  occurs. 

Comparison  of  Results.  First,  the  distribution  of  pressure 
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Figure  28  -  Convergence  Eeliavior  of  Potential  Function 
at  the  Midchord  as  Function  of  Free-Stream 
Mach  Number 
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coefficient  for  a  given  discretization  is  shown  in  Fig  29  as  a 
function  of  Mach  nunber.  The  distributional  shape  of  Cp  is  slightly 
different  than  predicted  from  linear  theory,  where  Cp  scales  with 
Mach  number  by  a  factor  (  |  -  •  The  departure  from  linear  theory 

is  not  significant  but  is  attributable  to  the  presence  of  the  nonlinear 
term  in  the  differential  equation,  which  begins  to  become  significant 
for  Mach  number  near  .  More  discussion  of  nonlinear  effects 

is  included  in  the  next  section.  Figure  30  compares  finite  element 
results  with  experimental  data  obtained  by  Knetchnel  (Ref  39)  for 
Me"  0.7c??  .  Good  agreement  exists  over  most  of  the  airfoil.  For 
comparison  purposes  a  finite  difference  solution  (Ref  40)  is  given  for 

-  0.7  (which  is  1%  lower  than  for  other  results).  In  addition, 
the  exact-linear,  thin-airfoil-theory  results  obtained  by  scaling  the 
incompressible  solution  are  shown.  The  farfield  boundary  locations  were 
the  same  for  both  the  finite  element  and  finite  difference  approxima¬ 
tions  (>max  =  =  1.5).  However,  the  number  of  nodes  used  for  the 

finite  element  method  was  more  than  an  order  of  magnitude  less  than  the 
number  of  grid  points  used  for  the  finite  difference  method  (i.e.  225 
nodes  compared  to  8514  grid  points ,  respectively) .  The  finite  dif¬ 
ference  solution  does  not  compare  with  the  experimental  data  as  well 
as  either  the  finite  element  or  the  exact  linear  solutions. 

Figure  31  compares  Cp  distributions  for  8  obtained 

from  the  present  finite  element  solution,  finite  difference  solutions 
(Ref  40),  and  from  linear,  thin-airfoil-theory  solutions.  The 
difference  in  pressure  distributions  between  the  solutions  from 
linear  theory  and  the  two  numerical  solutions  based  on  nonlinear  theory 
are  readily  detectable.  Differences  between  the  finite  element  and 


97 


Fijnire  29  -  Mach  Effects  on  C,,(y  -  0)  for 

a  6%-Thick  Parabolic -Arc  Airfoil 
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finite  difference  results  occur  for  the  entire  distribution.  These 
variations  are  due  to  the  inherent  difference  in  approximations  achieved 
by  each  method  and ,  perhaps  to  a  lesser  extent ,  the  finer  mesh  size 
used  for  the  finite  difference  solution. 

Transonic  Flow 

Iterative  Behavior  of  Solution  Scheme.  As  discussed  previously, 
the  iterative  solution  algorithm  does  not  converge,  as  it  does  for  sub¬ 
sonic  flow,  once  a  supersonic  bubble  larger  than  half  of  an  element 
forms  in  the  flowfield.  Figures  32  and  33  show  solution  results  after 
each  iteration  for  Mach  numbers  equal  to  0.84  and  0.86  respectively. 

The  solution  shown  in  Fig  32  converges  in  the  sense  of  eq  45  after  12 

iterations  for  a  tolerance  of  €  -  0.5*10  .  The  peak  Cp  is  slight- 

* 

ly  above  the  critical  value  of  Cp  -  0> 3Y7  >  but  the  supersonic  bubble 
only  engulfs  approximately  half  of  the  element  which  straddles  the 
point  of  maximum  thickness.  When  the  Mach  number  is  increased  to  0.86 
the  solution  fails  to  converge.  Figure  33  shows  what  happens  to  the 
solution  for  the  first  6  iterations.  The  solution  remains  symmetric 
(shockless)  and  after  3  iterations,  more  than  one  element  is  contained 
in  the  supersonic  bubble.  Solutions  for  further  iterations  (except 
iteration  8  which  goes  off  the  page)  are  shown  in  Fig  34.  The  spike  at 
the  midchord  continues  to  increase  until  iteration  9,  when  it  changes 
sign  and  creates  a  crevice  at  the  midchord.  Further  iterations  result 
in  the  fallen  spike  being  rebuilt,  which  eventually  occurs,  and  then 
the  process  is  repeated.  Thus,  the  solution  is  oscillatory,  and  although 
it  never  converges  neither  does  it  diverge.  If  the  Mach  number  is  in¬ 
creased  further,  violent  oscillations  occur  and  eventually  the  solution 
diverges . 
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Figure  32  -  Convergence  Behavior  of  Iterative  Figure  33  -  Divergence  Behavior  of  Itera 
Solution  Scheme  for  0.84  Solution  Schane  for  M  =  0.8< 


Figure  34  -  Further  Iterative  Behavior 
(Continued  from  Figure  33) 
for  1^,*  0.86 
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Figures  32-34  were  included  only  to  illustrate  the  behavior  of 
the  solution  at  the  on-set  of  transonic  flow  when  the  finite  elerrent 
solution  is  obtained  as  if  the  problem  was  elliptic.  When  the  co¬ 
efficient  of  ^>^y,  in  the  differential  equation  becomes  negative , 
the  differential  equation  changes  type  from  elliptic  to  hyperbolic 
which  forms  the  supersonic  bubble  in  the  flow.  Since  elliptic  dif¬ 
ferential  equations  are  fundamentally  different  than  hyperbolic  ones, 
then  formulations  and  algorithms  suitable  for  elliptic  equations  are 
not  expected  to  be  valid  for  hyperbolic  equations.  At  this  point  in 
the  solution  scheme,  the  "upwinding"  techniques  described  in  Chapter  IV 
were  enployed. 

Upwinding  Techniques .  The  upwinding  techniques  that  were  reported 
in  the  literature,  as  described  in  Chapter  IV,  were  tried  for  the  tran¬ 
sonic  flow  problem.  None  of  the  reported  techniques  were  able  to 
stabilize  the  iterative  solution  scheme,  and  a  discussion  of  the 
solution  behavior  for  each  of  these  techniques  is  omitted.  In  general, 
when  they  were  employed  one  of  two  things  happened.  First,  for  most 
cases  considered,  the  application  of  the  upwinding  technique  caused  the 
iterative  solution  scheme  to  diverge  imnediately.  Secondly,  for  a  few 
cases  where  the  solution  scheme  converged,  it  converged  to  a  solution 
that  was  not  physically  meaningful.  Due  to  the  inability  of  the 
reported  upwinding  techniques  to  stabilize  the  solution  scheme,  other 
upwinding  schemes  were  sought.  As  a  result,  the  new  proposed  up- 
winding  technique  was  developed. 

Recall  that  the  new  upwinding  technique  described  in  Chapter  IV 
modifies  the  finite  element  formulation  to  account  for  the  hyperbolic 
character  of  the  equation  for  those  elements  within  the  supersonic 
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bubble.  The  behavior  of  the  upwinding  procedure  is  governed  by  two 
parameters  R  and  V .  Equation  62  gives  the  modification  of  the 
finite  element  formulation  of  the  governing  equation  by  these  two 
parameters  for  elements  within  the  supersonic  or  hyperbolic  bubble. 

The  behavior  of  the  solution  as  a  function  of  these  two  parameters  is 
shown  in  Figs  35  to  39.  For  the  airfoil  considered,  the  numerical 
values  of  R  and  V  fall  somewhere  in  the  ranges :  0  —  &  —  0,  "i  and 
O  5  V  <  O.  £ 

The  behavior  of  the  solution  (pressure  distribution)  as  a  function 
of  upwinding  parameter  1/  for  fixed  values  of  M^and  parameter  R  is 
shown  in  Figs  35-37.  Although  a  different  value  of  R  or  rA^is  held 
fixed  in  each  of  these  figures ,  the  effect  of  If  on  the  solution  is  the 
same  in  each  case .  The  function  of  V  is  to  change  the  symmetry  of  the 
flow.  If  V  is  set  equal  to  zero  (corresponding  to  the  upwinding  scheme 
not  being  employed) ,  then  the  distribution  of  Cp  remains  syrimetric 
about  the  midchord,  as  shown  in  Fig  33.  Note  in  Figs  35-37  that  for 
increasing  values  of  V  the  distributions  are  skewed  downwind.  Pressure 
gradients  become  greater  on  the  downwind  side  of  the  peak  value,  in  the 
vicinity  of  an  expected,  weak  compression  shock.  Thus  parameter  \J 
not  only  alters  the  symmetry  of  the  pressure  distribution,  but  it  also 
appears  to  capture  the  behavior  of  weak  shocks  by  permitting  relatively 
large  discontinuities  of  velocity  to  occur  aft  of  the  peak  pressure. 
This  behavior  is  particularly  evident  in  Fig  37  for  1/  -  O^O. 

The  behavior  of  the  solution  (pressure  distribution)  as  a  function 
of  parameter  R  for  fixed  values  of  V  is  shown  in  Figs  38  and  39.  The 
distribution  of  pressure  coefficient  is  not  substantially  affected  by 
changes  in  parameter  R  except  near  the  peak.  In  this  region  increased 
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Figure  35  -  Effect  of  Upwinding  Parameter  U  for 
a  67c- Thick  Parabolic-Arc  Airfoil  at 
0.908  and  R  =  0 
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Figure  36  -  Effect  of  Unwinding  Parameter  U  for 
a  67-TViick  Parabolic -Arc  Airfoil  at 
M  =  0.908  and  R  «  0.30 
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Figure  37  -  Effect  of  Upwinding  Parameter  U  for 
a  -Thick  Parabolic-Arc  Airfoil  at 
=0.92  and  R  -  0.20 
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values  of  R  produce  larger  peak  pressures  which  means  aft  of  the  peak 
larger  gradients  of  pressure  are  produced.  Thus  parameter  R  also 
appears  to  be  contributing  to  the  numerical  mechanism  of  capturing  the 
weak  compression  shock.  Recall  that  the  function  of  R  is  to  retain 
part  of  the  nonlinear  term  on  the  left-hand-side  of  the  governing 
equation  and  to  subject  the  remaining  part  to  the  upwinding  operation. 
Judging  from  the  results  in  Figs  38  and  39  it  may  be  necessary  to  keep 
part  of  the  term  on  the  left-hand-side  of  the  equation,  which  means 
that  an  inversion  must  be  done  for  each  iterate.  If  the  entire  non¬ 
linear  term  is  treated  as  a  force  (i.e.  R  =  O  )  then  accuracy  may  be 
sacrificed.  The  value  of  R  for  the  cases  considered  range  between 
0.2  and  0.3  for  best  results. 

If  R  is  arbitrarily  selected,  say  ft  "  ft*  ,  for  given  values  of 
and  V,  then  one  of  three  things  could  happen.  First,  if  is 
picked  too  large  the  solution  will  either  diverge  or  oscillate  about 
some  solution  as  demonstrated  in  Figs  33  and  34.  For  this  situation 
not  enough  of  the  nonlinear  term  is  subjected  to  the  upwinding  operation 
to  permit  convergence  of  the  iterative  scheme.  Secondly,  if  R*  is 
picked  too  small,  then  too  much  of  the  nonlinear  term  is  altered  by 
the  upwinding  process.  As  a  consequence,  even  if  the  solution 
algorithm  converges,  the  solution  may  not  be  as  close  to  the  true 
solution  as  it  could  be  for  a  larger  value  of  .  Thirdly,  if  R&.is 
appropriately  chosen,  the  iterative  scheme  will  converge  to  an  accurate 
solution  which  correctly  describes  the  physics  of  the  problem  as 
shown  for  the  examples  in  the  next  section. 

In  general .  the  values  of  R  and  “IT  depend  upon  the  airfoil  shape , 
Mach  number,  and  the  angle  of  attack  for  lifting  airfoils.  When 
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values  of  R  and  V  are  correctly  chosen ,  then  convergence  occurs  rapidly ; 
usually  less  than  15  iterations  are  required.  If,  on  the  other  hand, 
they  are  not  correctly  chosen,  then  the  iterative  scheme  will  not 
converge.  After  a  couple  of  iterations  the  scheme  will  reach  a  point 
of  divergence  which  indicates  that  the  other  values  of  R  and  V  should 
be  selected  for  the  case  being  considered.  Unfortunately,  there  are 
no  known  analytical  expressions  which  select  the  best  values  of  R  and 
V  for  specific  airfoils  and  flow  cases.  The  values  can  be  determined 
iteratively  as  suggested  by  the  procedure  in  Chapter  IV.  This  proce¬ 
dure  first  requires  that  "ballpark"  values  of  R  and  V  be  found  which 
permit  convergence  of  the  iterative  scheme.  Next,  R  and  V are  "fine 
tuned"  to  select  the  "best"  possible  values. 

Comparison  of  Results .  Finite  element  solutions  for  transonic 
flow  over  a  6%  -  thick  parabolic-arc  airfoil  are  compared  with  experi¬ 
mental  data  obtained  by  Knetchel  (Ref  39) ,  finite  difference  calcu¬ 
lations  of  Olsen  and  Batill  (Ref  40)  and  also  those  of  Olsen  (Ref  60) , 
and  with  finite  element  solutions  of  Akay  (Ref  24).  For  -  0.9c8 , 

finite  element  solutions  for  the  pressure  coefficient  are  presented 
in  Figs  40  and  41  for  seven  discretizations  of  the  flow  domain  ftp,  as 
specified  in  Table  VI. 


Table  VI 

Domain  Discretization  Parameters, 

For  s  0.^0  8 


Grid 

NDX 

NDXA 

NDY 

Nodes 

1 

8 

5 

6 

98 

2 

8 

7 

7 

128 

3 

8 

7 

9 

160 

4 

6 

14 

9 

210 

5 

1 

6 

9 

13 

224 

6 

!  6 

i 

9 

1 

13 

224 

7 

6 

11 

i _ 

11 

216 

By  comparing  the  solutions  for  the  pressure  coefficients  given  in 
Fig  40  with  those  given  in  Fig  41,  a  significant  difference  in  behavior 
is  noticed  only  in  the  vicinity  of  the  expected,  weak  compression 
shock.  In  Fig  40  for  grids  1-3  at  a  location  eight-tenth  of  a  chord- 
length  frcm  the  leading  edge  (  =  0.8  )  ,  a  comparatively  large  jimp 

in  Cpi-i-c.)  occurs  between  adjacent  elements.  This  behavior  occurs 
naturally  in  the  solution  process  and  is  associated  with  the  occurrence 
of  a  weak  compression  shock  which  occurs  in  the  flow  domain  to  permit 
the  fluid  in  the  supersonic  bubble  to  return  to  subsonic  conditions. 

The  occurrence  of  the  jump  at  XL£  -0.8  was  not  forced  by  imposing 
any  shock  jump  conditions,  but  it  is  an  inherent  consequence  of  the 
potential  formulation  of  the  problem,  the  finite  element  approach 
selected  to  solve  the  problem,  and  the  use  of  the  new  upwinding  scheme 
governed  by  parameters  R  and  \S .  Since  the  velocity  potential  function 
is  being  solved  for,  and  since  elements  were  selected  which  insure  only 
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Figure  41  -  Finite  Element  Solutions, 
R  =  0.25  and  U  =0.50 
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continuity  of  the  velocity  potential  function  across  inter-element 
boundaries,  then  jumps  or  discontinuities  in  the  derivatives  of  the 
potential  function  (i.e.  velocities)  across  inter-element  boundaries 
are  inherently  expected  to  occur.  Fran  small-disturbance  theory,  the 
pressure  coefficient  is  proportional  to  the  velocity  in  the  x-direction, 
thus  a  corresponding  jump  of  pressure  between  elements.  The  signifi¬ 
cant  part  of  this  behavior  is  that  the  jump  is  of  sufficient  magnitude 
to  allow  the  flow  to  change  from  supersonic  to  subsonic  conditions  by 
crossing  an  element  boundary  line,  thus  an  ideal  shock  of  zero  thickness. 
This  particular  behavior  occurs  only  when  the  discretization  is  arranged 
in  such  a  manner  that  a  node  occurs  at  the  location  where  the  weak 
shock  would  like  to  impinge  upon  the  airfoil. 

The  solutions  presented  in  Fig  41  are  notably  different  in  the 
vicinity  of  xL£*o,e  For  these  solutions  the  discretization  does 
not  place  a  node  at  Xig  -  0.  8  ;  the  discretizations  straddle  that 
point  instead.  Since  the  shock  should  occur  at  Xl£  -  o.e  .  then  it 
must  occur  within  an  element  rather  than  on  its  boundary.  However, 
there  is  no  numerical  mechanism  for  the  shock  to  occur  within  an 
element,  as  there  is  between  elements.  Consequently,  two  jumps  are 
required  to  return  the  flow  from  supersonic  to  subsonic  conditions. 

For  these  discretizations  the  solutions  smear  the  shock  over  the  span 
of  one  element  and  two  successive  inter-element  boundary  locations. 

The  difference  between  solutions  given  in  Fig  40  and  those  given 
in  Fig  41  is  more  graphically  demonstrated  in  Fig  42,  where  solutions 
for  grids  3,  6,  and  7  are  compared.  The  symbol  at  the  peak  of  each 
distribution  indicates  the  location  of  the  expected  weak  shock,  which 
is  denoted  by  the  vertical  arrow  at  Xi>  £  ~  €>•  &  All  other  symbols 
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are  plotted  at  the  midpoint  of  the  interval  (i.e.  midway  between 
nodes)  as  described  previously.  Note,  as  the  elements  which  straddle 
the  shock  become  large  (i.e.  grid  6  as  compared  to  grid  7),  the  more 
the  solution  is  affected  in  the  vicinity  of  the  weak  shock.  It  appears 
that  smearing  the  shock  over  a  larger  element  tends  to  decrease  the 
magnitude  of  the  peak  pressure,  as  well  as,  affecting  the  magnitude  of 
the  jumps  required  to  return  the  flow  to  subsonic  conditions.  Although 
smearing  does  not  provide  the  best  approximation,  in  terms  of  describing 
the  physics  of  the  flow,  it  does  provide  an  acceptable  approximation 
without  the  use  or  need  of  shock  elements.  Better  approximations  can 
be  obtained  by  altering  or  "fine  tuning"  the  discretizations  to  provide 
solutions  similar  to  those  obtained  from  grid  3. 

The  new  upwind  method  appears  to  exhibit  characteristics  corrmon  to 
both  shock  "capturing"  and  "fitting"  techniques.  The  location  and 
strength  of  the  weak  shock  can  be  obtained  without  fitting  the  discreti¬ 
zation  in  the  sense  described.  However,  if  better  approximations  are 
required,  then  the  discretization  can  be  iterated  upon  to  fit  the 
element  boundaries  so  they  coincide  with  the  location  of  the  shock. 

From  the  approach  taken  in  this  study,  the  shape  of  the  shock  is 
restricted  to  a  vertical  line  segment  (i.e.  parallel  with  the  y-axis) . 

In  general,  the  shock  may  be  inclined  to  the  vertical,  although  it 
would  remain  straight.  Discretizations  could  be  devised  by  rotating 
segments  of  the  grid  to  fit  the  shock  between  adjacent  elements  for 
such  cases.  Neither  of  the  situations  represented  in  Fig  42  would  be 
possible  without  the  use  of  the  new  upwinding  technique.  Without  its 
use  the  pressure  distribution  would  remain  synmetric  with  respect  to 
the  midchord,  and  the  solution  algorithm  would  not  converge. 
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Convergence  is  possible  only  because  the  physics  of  the  flow  is  modelled 
by  excluding  the  downwind  influence  on  the  solution  at  upwind  nodes. 

Figure  43  compares  the  finite  element  solutions  for  no  smearing  of 
the  shock  with  the  finite  difference  solutions  and  with  experimental 
data.  The  two  finite  difference  solutions  correspond  to  a  constant 
grid  size  (Ref  40)  and  a  variable  grid  size  (Ref  60) ,  respectively. 

The  latter  solution  was  obtained  for  smaller  step  sizes  near  the  leading 
and  trailing  edges  and  in  the  vicinity  of  the  expected  shock.  The 
finite  element  solutions  were  also  obtained  from  both  constant  (grid  1) 
and  variable  element  sizes  (grids  2  and  3) .  For  the  variable  element 
discretizations,  smaller  (more  narrow)  elements  were  placed  near  the 
leading  and  trailing  edges.  In  addition,  solutions  were  attempted  for 
discretizations  that  were  refined  only  in  the  region  where  the  shock  was 
expected  to  occur.  This  discreti  stion  resulted  in  a  grid  where 
’’needle  like”  elements  were  placed  adjacent  to  elements  of  aspect  ratio 
near  unity.  The  solution  algorithm  for  such  discretizations  did  not 
converge  and  the  idea  of  using  variable  size  elements  to  assist  in 
’’capturing”  the  shock  was  abandoned.  It  is  believed  that  the  upwinding 
technique  is  applicable  only  for  elements  of  constant  size  throughout 
the  supersonic  bubble,  and  it  works  best  when  the  aspect  ratio  of  the 
elements  is  near  unity.  Any  substantial  departure  from  a  uniform  grid 
particularly  in  the  center  segment  (refer  to  Fig  13)  leads  to  conver¬ 
gence  problems  in  the  iterative  solution  process. 

Figure  44  shows  a  comparison  of  pressure  distributions  obtained 
from  two  finite  element  solutions  using  different  upwinding  techniques . 
The  finite  element  solution  presented  in  Fig  40  for  grid  3  is  compared 
with  the  solution  obtained  by  Akay  (Ref  24)  .  His  upwinding  technique 
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differed  from  that  used  in  this  study  as  described  in  Chapter  IV.  It 

z  viz. 

consisted  of  replacing  the  velocity  $eby  +  ,  where 

is  the  velocity  upstream  of  element  6  .  The  value  of  ©  was  taken 
to  be  0.20.  The  results  of  this  study  agree  with  the  trend  of  the  ex¬ 
perimental  data  and  the  finite  difference  solutions  much  better  than 
do  Akay's  predictions.  It  is  believed  that,  the  difference  in  upwind 
techniques,  rather  than  formulation  methods  is  what  accounts  for  the 
variation  between  the  two  Cp  distributions. 

Figure  45  compares  the  distribution  of  Cf(v=o)  with  experimental 
data  for  a  6%- thick  parabolic  arc  airfoil  at  ~  0.^2.  .  Domain 

Xlf  was  discretized  with  grid  7  described  in  Table  VI.  Upwinding 
parameters  of  R  “  0.2.0  and  V  -  0.  V0  were  required  to  obtain  conver¬ 
gence.  Essentially  the  same  behavior  observed  for  the  solution  for 
M  -  0  ^  applies  to  this  case.  The  shock  is  smeared  over  one 

CD 

element  and  the  pressure  distribution  agrees  with  the  trend  of  the 
experimental  data. 

In  general  the  new  upwinding  technique  gave  acceptable  results . 

A  number  of  other  upwinding  techniques  described  in  Chapter  IV  were 
tried.  None  of  these  gave  satisfactory  results,  although  each  tech¬ 
nique  did  alter  the  symmetry  in  the  flow.  In  most  cases  the  initial 
application  of  the  upwinding  technique  resulted  in  a  dramatic  change  in 
the  flowfield  velocities.  As  a  result,  the  solution  scheme  either 
diverged  after  one  or  two  iterative  steps,  or  it  oscillated  without 
regularity  about  some  solution.  The  possibility  exists  that  these 
upwinding  techniques  could  be  modified  to  provide  successful  application, 
but  no  such  alterations  were  found  in  this  study. 
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VII  Summary,  Conclusions ,  and  Recommendations 


The  Galerkin  Finite  Element  Method  was  used  to  obtain  approximate 
solutions  of  the  steady,  transonic,  small -disturbance,  velocity  po¬ 
tential  equation  for  flow  over  a  circular  cylinder  and  thin  airfoil. 

Circular  Cylinder 

Incompressible  Flow.  One  non-conforming  and  two  conforming  sector 
elements  were  used  to  solve  the  problem  of  incompressible  flow  without 
circulation  over  a  unit  circular  cylinder.  The  non-conforming  element 
was  rejected  because  the  solution  error  became  significant  as  element 
size  was  refined.  It  gave  accurate  approximations  only  for  coarse 
discretization  of  the  flowfield.  Velocity  distributions  calculated 
from  the  two  conforming  elements  agreed  well  with  the  exact  solution. 

In  general,  the  new  sector  element  gave  more  accurate  solutions  than 
the  other  sector  element ,  which  was  based  upon  a  conventional  bilinear 
polynomial  approximation.  The  new  element  was  developed  from  a  trial 
solution  using  rational  functions.  For  incompressible  flow  over  the 
cylinder  the  radial  size  parameter  of  the  new  element  did  not  have 
to  be  refined  to  achieve  convergence  as  element  size  approached 
zero.  The  only  refinement  required  for  convergence  was  the  angular 
width  of  the  element,  which  represents  an  improvement  over  the  conven¬ 
tional  conforming  element.  This  element  required  refinement  of  both 
size  parameters  with  the  aspect  ratio  held  fixed  (equal  to  one)  to 
achieve  convergence.  Consequently,  for  a  desired  degree -of -accuracy 
more  elements  and  larger  degrees -of- freedom  were  needed  for  the 
bilinear  element  than  for  the  new  element.  This  translated  into  more 
computer  core-storage  requirements  and  longer  processing  time  to  solve 
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a  given  problem. 

The  differences  between  the  two  conforming  elements  used  for 
this  problem  showed  the  advantage  of  using  elements  with  trial 
functions  which  resemble  the  expected  form  of  the  solution.  This 
idea  could  be  extended  to  the  airfoil  problem,  or  for  that  matter,  any 
problem,  providing  the  behavior  of  the  solution  is  known.  The  solution 
for  the  airfoil  problem  at  least  in  the  farfield  behaves  like  */i r  .  It 
may  be  possible  to  use  rational  trial  functions  (perhaps  raised  to 
some  power)  to  approximate  the  solution  and  achieve  an  element  with 
improved  characteristics . 

The  utility  of  the  new  element  is  not  fully  realized  for  solving 
circular  cylinder  problems,  but  can  be  realized  when  solving  airfoil 
problems.  The  sector  element  could  be  used  to  discretize  the  flow- 
field  fr an  one  chord-length  or  other  specified  radius  on  out  to  the 
farfield  boundary.  Since  the  trial  function  for  the  sector  element 
automatically  satisfies  the  infinity  condition  (as  r ->  eo  ) ,  then  the 
boundary  could  be  extended  as  far  as  necessary  to  approximately 
satisfy  the  boundary  condition  there.  It  is  also  possible  that  an 
infinite  element  could  be  developed  from  this  element  by  neglecting  the 
polynomial  terms  in  the  trial  function.  These  elements  could  prove 
quite  useful  in  this  respect  since  a  farfield  boundary  condition 
given  by  eq  50  would  not  have  to  be  used,  and  the  dreadful  nonlinear 
term  that  requires  integration  over  the  flowfield  would  not  have  to 
be  evaluated. 

Circulation .  The  two  conforming  sector  elements  were  used 
to  solve  the  problem  of  incompressible  flow  with  circulation. 

A  splitting  technique  was  employed  which  permitted  the  Kutta  or 
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stagnation  condition  to  be  enforced  in  a  convenient  manner  after 
the  elementary  or  component  solutions  were  obtained.  Both  elements 
predicted  results  accurately.  In  general,  the  new  element  was  more 
efficient  and  required  less  refinement  for  a  given  degree- of- accuracy. 
The  circulation  could  be  predicted  with  its  use  to  within  0.0037»  of 
the  exact  value  for  coarse  discretizations.  This  degree-of-accuracy 
is  possible  only  when  the  discretization  is  constructed  in  such  a 
manner  that  the  stagnation  point  is  located  at  the  angular  centroid 
of  the  stagnation  element.  As  the  stagnation  point  is  moved  toward 
the  edge  of  the  element  the  error  in  circulation  increases,  and 
reaches  a  maximum  when  the  stagnation  point  lies  on  a  node.  Velocity 
distributions  were  accurately  predicted  by  both  elements  whenever  the 
circulation  was  correctly  predicted. 

Compressible  Flow.  The  new  sector  element  was  used  to  solve  the 
small-disturbance  equation  for  compressible,  potential  flow  over  the 
cylinder.  The  equation  was  locally  linearized  by  an  iterative  solu¬ 
tion  scheme  which  converged  rapidly  for  subsonic  flows.  The  scheme 
failed  to  converge  for  transonic  flows  when  the  supersonic  zone 
engulfed  at  least  one  complete  element.  Predictions  of  the  critical 
Mach  number  and  subsonic  velocity  distributions  compared  well  with 
known  results. 

Thin  Airfoil 

Bilinear  rectangular  elements  were  used  exclusively  to  discretize 
the  flow  dcmain  for  the  airfoil  problem.  Cases  of  incompressible  flow 
and  compressible  subsonic  and  transonic  flows  were  considered  for  a 
nonlifting  symnetric  airfoil.  The  governing  nonlinear  equation  was 
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written  as  a  sequence  of  linear  equations  with  variable  coefficients 
and  was  solved  by  an  iterative  process.  The  iterative  solution  al¬ 
gorithm  converged  very  rapidly  for  low  subsonic  flows.  As  the  Mach 
number  approached  a  critical  value  the  rate  of  convergence  decreased. 
For  Mach  numbers  slightly  below  critical,  convergence  still  occurred 
in  less  than  twelve  iterations.  Convergence  continued  to  occur  as  the 
Mach  number  was  increased  until  more  than  half  of  one  element  became 
engulfed  in  the  supersonic  region.  At  this  point  the  iterative  scheme 
oscillated  about  some  solution  and  diverged  for  still  larger  Mach 
numbers . 

To  achieve  convergence  for  transonic  flows  (mixed  elliptic- 
hyperbolic)  ,  the  finite  element  formulation  had  to  be  modified  to 
account  for  the  proper  zones  of  influence  for  those  elements  within 
the  supersonic  (hyperbolic)  region.  Several  modifying  methods  in¬ 
cluding  those  reported  by  other  investigators  were  tried.  All  but 
one  method  failed  to  produce  the  desired  result.  This  method  is  a 
new  upwinding  technique  governed  by  two  parameters  which  exclude  the 
influence  of  iterative  downwind  forces  on  the  solution  at  upwind 
nodes.  The  new  upwinding  technique  not  only  keeps  the  iterative 
solution  scheme  from  diverging,  but  it  also  captures  the  weak  com¬ 
pression  shock  which  forms  in  the  flowfield. 

Accurate  approximations  for  pressure  distributions  were  obtained 
for  all  flow  regimes  from  incompressible  to  transonic  flows.  The 
simplifying  assumptions  and  approximations  used  in  this  study 
represents  an  exceedingly  basic  approach  to  an  extremely  complex 
nonlinear  problem.  One  of  the  original  intents  was  to  demonstrate 
whether  a  simple  approach  was  acceptable.  The  finite  element 
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techniques  used  in  this  study  represent  approaches  that  are  about 
as  basic  as  one  could  take.  Excluding  the  use  of  linear  triangular 
elements,  it  would  be  difficult  to  conceive  of  a  more  fundamental  way 
to  solve  the  nonlinear  small-disturbance  equation  than  the  way  it  was 
done  in  this  study.  More  complicated  problem  formulations  as  well  as 
higher  orders  of  approximation  could  be  used.  Some  of  these  approaches 
are  described  in  the  remaining  section. 


Reconmenda t  ions  for  Future  Work 

Further  development  of  finite  element  methods  for  airfoil  analysis 
could  take  numerous  directions.  Several  extensions  to  the  present 
effort  could  be  made,  as  well  as,  additional  approaches  to  the  problem. 
Future  work  will  undoubtedly  be  conducted  in  the  areas  of:  problem 
formulation,  element  development,  application  of  existing  higher-order 
elements,  singular  treatment  of  the  leading  edge,  unsteady  analysis,  and 
special  treatment  for  the  mixed  or  transonic  problem.  Each  of  these 
areas  will  be  discussed  briefly  with  specific  suggestions  for  develop¬ 
ment  programs. 

Problem  Formulations .  For  steady  flow  three  additional  inviscid 
formulation  techniques  could  be  investigated.  They  are  the  velocity 
formulation  for  small-disturbance  theory  and  both  velocity  and  potential 
formulations  for  large  disturbances. 

The  small-disturbance  potential  equation  used  in  this  study  can 
alternatively  be  written  as  two  equations  in  terms  of  the  disturbance 
velocities  as 

ut%  +  -  O  (66) 

and 

-  y*  =  O  (67) 

Boundary  conditions  are 

/  \  cJ"f 

(' +  u )  77  (68) 
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and 


(u,w  )  — >  O  as  Y"  — (69) 

These  conditions  are  Dirichlet  conditions  which  are  imposed  upon  the 
trial  functions.  The  manner  in  which  the  boundary  conditions  are 
treated  using  this  formulation  technique  differs  from  the  small- 
disturbance  potential  formulation  used  in  this  study.  The  present 
boundary  conditions  are  of  the  Neumann  type  which  are  not  imposed  on 
the  trial  function.  They  are,  instead,  enforced  indirectly  by  the 
formulation  procedure. 

The  velocity  formulation  with  the  use  of  linear  elements  as 
used  in  this  study  would  provide  pressure  distributions  along  the  air¬ 
foil  contour  that  are  linear  and  continuous,  without  jumps  between 
elements.  This  would  be  an  improvement  over  the  potential  formulation 
which  results  in  pressure  distributions  that  are  step-functions  for 
linear  elements  and  step-linear-functions  for  quadratic  elements.  To 
achieve  continuous  pressure  distributions  from  the  potential  formu¬ 
lation  Hermite  polynomials  would  have  to  be  used  as  approximation 
functions.  Along  with  the  improved  accuracy  of  the  velocity  formula¬ 
tion  over  the  potential  formulation  comes  the  disadvantage  that  twice 
as  many  degrees -of- freedom  are  needed  for  linear  element  discretiza¬ 
tions.  Two  unknowns  U-  and  \T  are  required  as  nodal  parameters 
instead  of  the  single  parameter  . 

The  other  problem  formulations  which  should  be  further  investi¬ 
gated  are  the  full-potential  and  alternatively  full-velocity  formula¬ 
tions.  The  full-potential  equation  for  inviscid  compressible  flow  is 

130 


t 


given  by  eq  15.  This  equation  could  be  cast  into  velocities  U  and 
\T  and  used  in  conjunction  with  the  irrotational  condition  These 
formulations  would  not  be  limited  to  slender  aerodynamic  bodies,  since 
large  disturbances  are  not  excluded  by  assumption.  However,  non¬ 
linearity  is  more  severe  than  for  the  small-disturbance  equation  which 
will  create  numerical  difficulties.  As  a  first  approximation,  at 
least  for  subsonic  flow,  the  equation  could  be  cast  into  an  iterative 
Poisson  equation.  Perhaps  an  upwinding  technique  similar  to  that  used 
in  this  study  could  also  be  tried  for  supercritical  flow.  Other 
solution  possibilities  exist  which  should  be  investigated. 

Higher  Order  Elements.  A  direct  extension  of  the  present  study 
is  possible  by  using  higher  order  elements  to  approximate  the  solution. 
Initially,  only  those  elements  along  the  airfoil  contour  should  be 
made  of  higher  order  with  all  other  elements  remaining  bilinear.  This 
arrangement  would  permit  an  improvement  in  the  description  of  the 
pressure  distribution  along  the  airfoil  contour  without  significantly 
increasing  the  total  number  of  nodal  parameters.  As  a  further  extension, 
all  elements  could  be  of  higher  order.  Results  from  these  solutions 
could  be  compared  with  those  of  this  study  to  determine  the  rate  of 
convergence,  (p  -type  convergence) . 

Boundary  Conditions.  Extensions  to  the  present  solutions  could  be 
made  by  applying  the  tangential  boundary  condition  along  the  airfoil 
surface  rather  than  along  as  done  in  thin-airfoil  theory. 

This  improvement  is  not  expected  to  alter  results  very  much  for 
thin-airfoils ,  but  would  be  necessary  for  thick  ones.  A  nunber  of 
approximations  are  possible.  First,  the  actual  boundary  shape  could 
be  approximated  by  suitably  shaped  elements  with  the  boundary  condition 


applied  along  the  appropriate  element  edge.  Secondly,  isoparametric 
elements  with  one  side  shaped  to  approximate  the  contour  shape  are 
possible.  Of  these  types  several  possibilities  exist  depending  on 
the  order  of  the  element  chosen  and  the  actual  contour  shape. 

Arbitrary  Thin-Airfoil  Shapes .  For  the  present  approach  a  com¬ 
putational  subroutine  needs  to  be  developed  to  extend  the  present 
analysis  for  an  arbitrarily  shaped  airfoil.  Two  extensions  are  re¬ 
quired  to  accomplish  this  goal .  First ,  the  nonlinear  form  of  the 
governing  differential  equation  will  have  to  be  included  for  lifting 
flow.  This  addition  will  also  require  the  inclusion  of  an  upwind ing 
technique  when  the  flow  is  transonic  These  additions  will  follow  the 
analysis  presented  for  the  syrrmetric  transonic  case  identically,  ex¬ 
cept  the  lower  half-space  will  be  included.  As  a  consequence,  inte¬ 
gration  along  the  lower  airfoil  surface  will  have  to  be  included. 

Secondly,  an  algorithm  will  have  to  be  added  to  the  subroutine  to 
compute  the  finite  element  vectors  and  matrices  that  come  from  the 
boundary  integrals  along  the  airfoil  profile.  These  computations  will 
have  to  be  made  from  normalized  airfoil  profile  coordinates.  Since  the 
boundary  integrals  depend  on  the  airfoil  slope  then  a  procedure  to 
accurately  describe  the  slope  from  profile  coordinates  is  needed.  For 
a  family  of  airfoils  which  can  be  defined  by  an  equation  with  perhaps 
variable  coefficients  the  integral  can  be  computed  exactly  (numerical 
integration  may  be  easier).  A  separate  routine  would  be  needed  for 
each  equation  (family)  type  with  the  coefficients  as  input  parameters. 

New  Elements.  An  interesting  study  would  be  to  attempt  develop¬ 
ment  of  new  elements  for  compressible  flow  problems  which  have 
properties  similar  to  the  new  element  used  in  this  studv  for 
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incompressible  flow  about  a  circular  cylinder.  The  first  attempt  might 
be  for  the  circular  cylinder;  although,  the  thin-airfoil  may  be  easier, 
due  to  the  decrease  in  nonlinearity.  Trial  functions  should  be  at¬ 


tempted  which  resemble  the  form  of  the  solution.  For  the  cylinder 

% 

problem,  forms  with  CM  might  be  attempted  where  4  is  variable.  For 

%  % 

the  airfoil  problem,  forms  with  and  could  be  tried.  Per¬ 

haps,  terms  which  vary  as  the  farfield  expression  [i.e. 


could  also  be  included. 


Shock  Elements .  An  area  which  has  not  been  successfully  treated 
with  finite  element  methods  for  2-D  problems  is  the  correct  treatment 
of  shock  waves.  Shock  fitting  techniques  could  be  developed,  but  a 
better  and  more  elegant  technique  would  be  to  develop  a  ’’special"  shock 
element  to  capture  both  the  location  and  strength  of  the  shock.  Such 
an  element  would  have  to  be  constructed  from  discontinuous  functions. 


Since  the  location  of  the  shock  is  not  known  a  priori,  the  element 
would  have  to  be  sufficiently  general  to  exclude  shocks  if  they  do  not 
appear  and  describe  their  strength  and  relative  location  within  the 
element  when  they  do  appear.  These  "special"  elements  would  be  used 
only  in  the  region  of  anticipated  shock  occurrence. 

Singularities .  For  airfoils  with  sharp  leading  edges  or  other 
slope  discontinuities  along  the  contour  there  will  be  local  singularities. 
Rather  than  discretizing  these  regions  with  extremely  small  elements  to 
describe  the  rapid  change  of  behavior,  a  local  subregion  could  be 
isolated  to  treat  the  singularity  separately.  The  use  of  "special" 
shape  functions  which  describe  the  singularity  near  the  singular  point 
could  be  included  in  the  finite  element  analysis.  This  technique  has 
been  done  in  applications  of  other  singular  problems,  such  as  the  stress 
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concentration  near  a  crack  tip. 

Unsteady  Analysis.  Chan's  et.  al.  work  is  the  only  extensive  un¬ 
steady  analysis  to  appear  in  the  literature.  Additional  finite  element 
solutions  for  this  problem  are  required.  A  number  of  approaches  and 
formulation  procedures  are  possible.  One  particular  approach  which 
should  be  tried  for  the  small-disturbance  equation  is  to  assume  that 
the  potential  can  be  expressed  as 

t)  =  (70) 


The  spacial  shape  functions  are  the  same  as  those  used  for 

the  steady  problem.  What  is  different  is  the  nodal  parameters  4.  (t) 
are  functions  of  time.  Instead  of  obtaining  a  set  of  algebraic 
equations  to  solve,  a  set  of  ordinary  second-order  differential  equations 
will  result  for  the  nodal  parameters.  For  small  oscillations  about  the 
mean  steady  position  these  equations  will  be  linear  and  can  be  solved 
by  existing  methods. 
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Appendix  A 


Shape  Functions  and  Elemental  Equations 
for  Sector  Element  (1) 

This  appendix  contains  a  description  of  how  shape  or  basis 
functions  are  derived  from  assumed  trial  functions.  The  description 
primarily  pertains  to  elements  used  to  solve  the  potential  flow  pro¬ 
blems  formulated  in  Chapter  III  for  the  circular  cylinder.  This  appendix 
also  contains  the  derivation  of  finite  element  equations  for  the  symme¬ 
tric  flow  case  when  sector  element  (1)  is  used. 

Shape  Functions 

In  general,  when  a  continuum  problem  is  solved  by  the  Finite  Element 
Method,  the  continuum  is  divided  into  a  finite  number  of  elements  which 
are  connected  at  discrete  points  situated  on  their  boundaries  called 
nodes.  The  continuum  which  lias  infinite  degrees-of-freedom  is  replaced 
by  a  finite  number  of  unknown  nodal  parameters.  These  parameters  are 
nodal -point  values  of  the  solution  function  or  its  derivatives,  de¬ 
pending  upon  the  complexity  of  the  approximation  desired.  Between  nodes 
(i.e.  inside  an  element)  the  solution  function  is  approximated  by  an 
assumed  functional  relationship  which  can  be  expressed  in  terms  of  the 
unknown  nodal  parameters.  The  approximating  function  (trial  function) 
can  be  written  as  a  linear  combination  of  the  shape  functions  and  the 
unknown  nodal  parameters.  For  example,  consider  the  incompressible  flow 
problem  described  in  Chapter  III.  An  assumed  trial  function  for  the 
sector  element  shown  in  Fig  3  can  be  written  as 


where  Nj  (n$)  are  the  basis  functions  and  the  nodal  value  of  the 
solution  at  node  J  .  The  repeated  index  J  implies  summation  fran 
j-  It  is  immediately  apparent  that  the  shape  function  should 


have  the  property  that 


(A- 2) 


where  £'•  is  the  Kronecker  delta.  For  simple  elements  it  is  often 

J 

easy  to  select  a  functional  form  for  the  solution  and  witn  the  use  of 
the  above  property  simply  write  down  the  basis  functions .  To  insure 
that  the  guessed  form  produces  a  conforming  element,  the  required  con¬ 
tinuity  conditions  have  to  be  satisfied.  For  more  complicated  approxi¬ 
mations  the  guess  technique  may  not  work  very  well,  and  a  more  direct 
method  will  have  to  be  used  (Ref  59).  When  using  the  direct  method,  the 
assumed  trial  function  must  be  written  with  as  many  undetermined  con¬ 
stants  as  there  are  nodal  parameters  in  the  element.  It  could  be  ex¬ 
pressed  as 


<t>(r,e)=  G 


(A- 3) 


where  (f.»)  are  independent  variables  and  the  Cis  are  undetermined  con¬ 
stants.  What  is  done  for  the  Finite  Element  Method  is  to  express  these 
constants  as  functions  of  the  element  geometry  and  the  unknown  nodal 
parameters  4;Aj-  This  process  is  started  by  first 

writing  eq  A- 3  at  each  node  to  give  a  system  of  n  equations  expressed  as 


<t>i  -  Gij  aj 


(A-4) 


The  constants  are  expressed  in  terms  of  the  nodal  parameters  by 


4 


inverting  Gt'j  to  give 


(A-5; 


Substituting  the  constants  into  the  trial  function,  eq  A- 3,  and  col¬ 
lecting  terms  will  produce  an  expression  in  the  form  of  eq  A-l. 

To  illustrate  the  direct  process,  consider  trial  function  (1)  for 
the  sector  element  shown  in  Fig  3.  Trial  function  (1)  is  given  in 
Chapter  III  by  eq  35  as 


(f>(r,e)  =  a,  +  b, 


(35 


Expressing  at  each  of  the  nodes  gives  a  system  of  4  equations  for  the 
4  unknown  nodal  values  of  (f>  .  These  equations  are  expressed  in  matrix 
form  as 
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Inverting  this  matrix,  substituting  the  constant  into  eq  35,  and  col¬ 
lecting  like  terms  of  (j) •  will  produce  the  desired  form.  The  resulting 

•J 

shape  functions  are 
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where  and  *=  Vft  as  defined  in  Fig  3.  It  is  easily 

demonstrated  that  the  shape  functions  satisfy  (Vj } 

However ,  further  examination  reveals  that  continuity  of  across  inter- 
element  boundaries  is  guaranteed  only  along  radial  lines.  Along  cir¬ 
cumferential  lines  ^  is  continuous  only  at  the  nodal  points.  The  use 
of  this  element  for  the  second  order  problem  of  incompressible  flow 
over  a  circular  cylinder  as  formulated  in  Chapter  III  will  give  a  non- 
conforming  approximation. 

Elemental  Equations 

The  elemental  equations  for  incompressible  flow  over  a  circular 

e  «  e 

cylinder  are  expressed  by  eq  30  as  f°r  (c,;=  I,-";  v). 

£  J  J  J 

Stiffness  matrix  given  by  eq  31  as 

\) 

Kij  =  Nj,r  +  r 1  y  rdrde  Oi 

£ 

0 

Each  element  of  K*  *  is  calculated  by  substituting  the  shape  functions 

<] 

p 

and  performing  the  required  integration.  For  example,  the  K,,  com- 

e 

ponent  of  K{j  is  obtained  by  setting  t ;  j  s  /  ,  and  substituting 


into  eq  31. 

The  result  is 
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where 


Integrating  and  substituting  the  limits  gives 
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Note  that  ku  depends  only  upon  the  two  elemental  geometric  parameters 
e(  and 

& 

The  remaining  elements  of  KXare  obtained  in  a  similar  manner  to 
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symmetric  stiffness  matrix  of  the 

form 
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The  constants  CL,  b,  C  are  given  by 


d  = 

PC**')  1 

2(o<-l) 
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9  (*<-•> 
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jBv 
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The  force  vector  is  given  by  eq  39  as 

-f,  -  N*  Co5@J 


r 


Substituting  the  shape  functions  eq  A- 7  into  this  expression  and  inte¬ 
grating  gives 


146 


4 


Jk. 
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The  constants  A^C.D  are  given  by 
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Velocity  Distribution 

The  dimensionless  tangential  perturbation  velocity  U  ©0,0) 
along  the  contour  of  the  cylinder  is  calculated  from  the  potential 

e 

function.  In  element  £,  ,  is  given  by  eq  28 

4>V,0)  =  Nt*  (r,  e)  4>t-  (28) 

e 

Velocity  IA  in  each  element  is  defined  from  the  potential  function  by 

o 
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U 
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Along  the  cylinder  surface  fr=i) 

u«0,e)= 


(A-12) 


Substituting  the  shape  functions  eq  A-7  into  eq  A-12  and  evaluating 
at  I  gives 
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Appendix  B 


Shape  Functions  and  Elemental  Equations  for 
Sector  Element  (2) 


Shape  Functions 

Consider  the  sector  element  shown  in  Fig  3  and  trial  function  (2) 
given  by 

e 

-  0lz  4-  bzr  +  C2  ©  +  06) 


From  Ref  (27)  this  approximation  can  be  written  in  the  form 

j  c  / }  }  4  where 
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The  shape  functions  could  also  be  obtained  by  the  direct  procedure 
described  in  Appendix  A.  It  is  easily  demonstrated  that  the  shape 
functions  provide  continuity  of  4  across  inter-element  boundaries .  Thus , 
element  (2)  is  a  conforming  element  for  the  finite  element  formulations 
described  in  Chapter  III. 
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Elemental  Equations 

The  elemental  equations  for  incompressible  flow  over  a  circular 
cylinder  are  expressed  by  eq  30  as  (p .  tc*  ^ j  ^  ^  '  *  fj  4*)  * 

Stiffness  matrix  is  given  by  eq  31  as 

e  6»  fy> 

Kq  =  J  J  (K>  Njf„  +  71  Nse  qj*^)rclrc)0  (31' 

0a 

e 

Each  element  of  Kp  is  calculated  by  substituting  the  corresponding 

)  e  e 

shape  functions  and  integrating.  For  example,  the  KM  component  of 

%  f 

is  obtained  by  setting  t"  j  -  |  and  substituting  eq  B-la  into  eq  31 
The  result  is 


r. 


(B-2) 


(B-3) 


Integrating  and  substituting  the  limits  gives 

«  i  -v<__  + 

K"  '  6<y- o  0 

which  depends  only  upon  the  two  elemental  geometric  parameters  o<  and 


The  remaining  elements  of  Kij'  are  obtained  in  a  similar  manner  to 
give  a  symmetric  stiffness  matrix  of  the  form 


Kb  = 
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where 


K{  z  a  +  «zc  -  e 
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constants  £  •  *  ^  £  in  the  above  expressions  are  given  by 
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The  force  vector  +•  is  given  by  eq  39  as 
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Substituting  the  shape  functions,  eq  B-l, into  this  expression  and  inte¬ 


grating  gives 
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1  -  11  it  y 

e 

The  above  expressions  define  -Pj  for  each  element  &  which  borders  the 
cylinder  surface.  For  all  other  elements  +7  is  identically  zero. 

Velocity  Distribution 

The  dimensionless  tangential  perturbation  velocity  U&  (  I  ; 
along  the  contour  of  the  cylinder  is  calculated  from  the  potential 
function.  In  element  is  given  by  eq  28 

(p\ a)  -  M;  (  0  ©)  4?  (28) 

e 

Velocity  in  each  element  is  defined  from  the  potential  function  by 

(4(r,e)=  N;i&(y-,b)4k 

Along  the  cylinder  surface  (Y"  —  I  ) 

MeO, e)  -  m,> 0,^)4  <B-» 

For  each  of  the  conformal  elements  this  expression  reduces  to 

UeO,e)  -  “  $  J  (b-7) 

Within  element  Q.  the  tangential  velocity  along  the  contour  is  a 
constant.  Thus,  between  elements  there  will  be  a  jump  in  L(e .  This 
implies  the  velocity  distribution  is  a  step  function  along  the  contour. 
The  tangential  velocity  defined  by  eq  B-7  for  element  (2)  is  also 
the  tangential  velocity  for  element  (3)  described  in  the  next  appendix. 
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Appendix  C 


Shape  Function  and  Elemental  Equations 
for  Sector  Element  (3) 


Shape  Functions 

Consider  the  sector  element  shown  in  Fig  3  and  trial  function  (3) 
given  by 

4>e{y,&)  -  0.3  +  +  07) 


This  approximation  can  be  written  in 
where 


e  e 

the  form  <j)  -  Nj  for  (  j  =  l; 


y<k/r  +  i /< 


N,  = 

i-  'K 

(C-la) 

fair  -l 
•  H  - 1 

©b-  e 

(C-lb) 

N3  c 

i/*-l 

#  -  ©«_ 

(C-lc) 

r 

fft./v’  -  '/*_ 

(C-ld) 

'I* 

These  shape  functions  are  determined  from  the  trial  function  by  the 
method  described  in  Appendix  A.  They  can  also  be  deduced  from  the 
shape  functions  in  Appendix  B  by  inverting  and  .  This  element 
represents  a  new  element  which  to  the  author's  knowledge  has  not  been 
used  to  solve  an  aerodynamics  problem. 

Elemental  Equations 

Incompressible  Flow.  The  elemental  equations  for  incompressible 


flow  over  a  circular  cylinder  are  expressed  by  eq  30  as 


K--  £  =  fj  (4j  =  ■  •'  v; 


(3C 


Stiffness  matrix  )('•  is  given  by  eq  31  as 

J 


K*:  =  J  +  VWrcJ 


9 


(31 


Each  element  of  K  is  calculated  by  substituting  the  corresponding 

3  e  e 


shape  functions  and  integrating.  For  example  the  K|,  component  of  K ,'y 
is  obtained  by  setting  t'sj  s  {  and  substituting  eq  C-2a  into  eq  31 
The  result  is 


Integrating  and  substituting  the  limits  gives 


a/  -  3 


-h 


A  <* 


K"  "  U°<'i)  sc*- o' 


(C-2 


which  depends  only  upon  the  two  elemental  geometric  parameters  <X  and 

»■ 

The  remaining  elements  of  ftj  are  obtained  in  a  similar  manner  to 
give  a  synmetric  stiffness  matrix  of  the  form 


The  force  vector  -pt-  is  given  by  eq  39  as 


ew 


f .  “  |  Nt-  toSB  cJS 


(39) 


f-  rv=  i 


Substituting  the  shape  functions,  eqs  C-l,  into  this  expression  and  inte¬ 
grating  gives 

cose*  ^ 

? 

o 


f?  = 


-  Sill  - 


Slrt  4- 


tosBb~  CoS  &A 


(C-4) 


The  above  expression  defines  for  each  element  e  which  borders  the 

n,® 

cylinder  surface.  For  all  other  elements  -ft'  is  identically  zero. 

Compressible  Flow.  For  compressible  flow  over  the  cylinder  the 
elemental  equations  are  given  by  eq  44  as 


[Kcj  +  ®  fc  +  ?.■<<« 


(C-5) 


Matrix  |^*»  and  vector  £t*  are  the  same  as  defined  for  the  incompressible 
problem .  Matrix  U:(r)  comes  from  the  nonlinear  term  in  the  governing 
differential  equation,  and  is  given  by  eqs  47  and  48.  It  can  be 
expressed  as 

L.-jC f)  '  -*Uj  +  K”0+*)](  A;;* a.';  -Cq  tP.y) 

(C-6) 
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A;;  '  ff  c°*e  'V'V  rd<-de 
J  H-e 

B;j  =  n.>  Wj,e  rd.'ds 

Sie 

Cl\  -  (*,<  Nj,e  +  N.>  NjV)  r  J  *d  • 

r  .Sin  e  CQ-S  ^  |  . 

Oti  =  J  7 Nt^V  de 


r=r«.=  i 


Substituting  the  shape  functions,  eq  C-2,  into  each  of  these  expressions 
and  integrating  gives 


ot  +  l 


a,  -  a, 

a, 


^2.  “ 

-  A2 

h3 


where 


(C-7) 


a,-[(el-2©fce-f6t-v2)^i^  t -f- 


-  ■  ■  mmfr 


-« J»> 


n* 


i 


<!<->) 


1  ~3«<  .  -jz  a 


Matrix 


becomes 


where 


i 


i 


* 


U-^4 


a-  ©t> 


cozzs  - 


Sir\  2& 


vrznvirrjr.  \  i-v EE  g  S3S 


c,  r 


-  j_ 


f  ^b+6a  _  €?b  +© 

L  s  8 


6a  -z& 


cos  z a 


©, 


31 


*  20,1  | 

fi  J  l 


CM  :  -^[-^  +  C0SZ©  -  \ 


The  matrix  Dj ;  comes  from  the  boundary  term  which  becomes 

J 


°‘j  = 


-c,  O 
o  o 
o  o 


o 

o 

o 

o 


o 

o 

cL 


(C-10) 


The  vector  ^  ((f)*)  also  ccmes  from  the  nonlinear  term  in  the  governing 
differential  equation  and  is  given  by  eq  49  as 

£ 

3,( P)  -  m^[|  +  k" (!+■»)]  \  Cose  N;  |  ae 

•r-^.= » 

Substituting  the  shape  functions,  eq  C-2,  into  eq  49  and  integrating 
gives 


V 


3;(0-  *£,[*  +  I 


o 

o 


(C-ll) 


5v 

L  J 
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The  tangential  perturbation  velocity  UqC  1,6)  is  derived  in 
Appendix  B  for  element  (2) ,  and  is  given  by  eq  B-7 .  For  element  (3) 
the  velocity  is  given  by  the  same  equation. 


Appendix  D 


Finite  Element  Equations  for  Flow  Over 
an  Airfoil  for  a  Bilinear,  Rectangular  Element 


Shape  Functions 

The  shape  functions  for  the  bilinear  rectangular  element  shown  in 
Fig  46  are  given  by 


N- =  -$-(>  +  (;  ()('  + 1-' *>)  <D-1) 


where  i  —  I ;  #  #  *  Y 

nates  of  node  c  . 


Coordinates 


are  the  local  nodal  coordi- 


Elemental  Equations 

The  weak  solution  of  the  governing  differential  equation  written 
in  elemental  form  is  expressed  by  eq  58  as 

[(i-Mi)  a.-j  e.j  -  m*  Csj(^) 

+  rvt  Dq  +  mI  42  ®-2) 


where 


I 


3 

\ 


I 


L 

i 

r 

r 


y 


B<j  =  f(  Nj,v  d*d* 

jsi.  e 

CiiW  =  ((  NX.<  N.'.x  NjJlt  Jx  Jv 


D‘i  dx  N‘  ^j’<*  I 

sn*  v=o 


M<  N;,*|  J* 

y- O 


4x  M>  I d 


x 


The  local  coordinates  <*,«!>  are  related  to  the  global  coordinates 
C*^)  by  the  transformation  equations 


X  ~  )^c 
a 


^  - 


i_-Jc 

b 


(D-l) 


Using  these  expressions  the  elemental  equations  can  be  written  in 
local  coordinates  to  give 


»  \ 

-i  - 1 
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N;,,  <J*  J«i 


J  Nk.*  N*  Nj,j| 

*1  --I 

I 

fc» 

-  l 

^*“1 


The  first  three  matrices  are  dependent  only  on  the  size  of  the  element 
and  are  not  related  directly  to  any  particular  airfoil  contour.  Sub¬ 
stituting  the  shape  functions,  eq  D-l.  into  these  expressions  and  letting 
t  -  i*  -  |  will  produce  the  following  results: 


- i 


1z  -  7z 

I  -  / 


Matrices  Dt‘j  and  £t’j  and  vector  *Pt'  are  directly  dependent  upon 
the  airfoil  contour  being  considered.  These  quantities  are  evaluated 
only  for  elements  that  share  a  corrrnon  boundary  with  the  contour.  Thus 
for  most  elements  in  the  domain  these  qualities  are  zero. 

Symmetric  Flow.  Consider  a  parabolic-arc  airfoil  with  thickness 
distribution  given  by 

•j  =  fcx)  =  r(  i-  vx1)  ;  jx 1 1  ‘It.  <p-v 

For  this  airfoil  the  remaining  quantities  in  eq  D-2  become 


o  o  o 

o  o  o 

o  a  -  a 

o  f 


o 


8T 


c 

o 

o 


o 

o 

o 

o 


O  D 
O  O 

H  -  H 


167 


Pressure  Distribution 


For  steady  small  disturbance  theory  the  coefficient  of  pressure 
is  given  by 

Cp ~  -  2.  u 


The  velocity  in  element  £  is  calculated  from  the  assumed  solution  for 
the  potential  given  by  eq  57 


(57) 


From  thin  airfoil  theory  the  pressure  coefficient  is  evaluated  along 
c  O*  .  Thus  for  element  £  which  borders  the  airfoil  contour  the 
elemental  pressure  coefficient  becomes 


o+)  =  -2  (X/y=o'fJ 


Using  the  transformation  equations  between  the  (X.V)  and 
coordinates,  given  by  eq  D-3,  and  substituting  eq  57  gives 


C  t,i) 


Cp(x.*8o*)*  ft(  £-<£) 


Within  element  Q  the  pressure  is  a  constant  value  which  means  there 
are  "jumps"  in  pressure  between  elements  along  the  airfoil  surface.  This 
implies  that  the  pressure  distribution  along  ^  :  o  +  is  a  step  function. 
For  most  of  the  distributions  of  the  pressure  coefficient  given  in 
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Chapter  VI  the  value  of  pressure  coefficient  computed  from  eq  D-5  is 
plotted  at  the  midpoint  of  the  element. 

Mach  Number 

The  local  Mach  number  is  defined  by  V/flLj  where  is  the  local 
speed  of  sound  given  by 


The  Mach  number  becomes 


The  velocity  V  is  given  by  +  elernent  £ 

( y/vcof  =  0+  +  l&'jf 

Transforming  to  local  coordinates  (5,^)  gives 

(v/u =  [  I  +  va  { tf)  +  0-  -£)}]*' 

(>- ati- $)}f 

Substituting  this  expression  into  eq  D-6  will  define  the  Mach  number 
in  element  &  .  For  the  bilinear  element  Mach  number  will  be  discon¬ 
tinuous  across  inter-element  boundaries. 
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